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SUMMARY 


A finite -difference procedure for computing the. turbulent, swirling, 
compressible flow in axisymmetric ducts is described. Arbitrary distributions of 
heat and mass transfer at the boundaries can be treated and the effects of struts, 
inlet guide vanes, and flow straightening vanes can be calculated. The calculation 
procedure is programmed in F0RTRAN IV and has operated successfully on the UNIVAC 1108, 
IBM 360, and CDC 66 OO computers. The analysis which forms the basis of the procedure, 
a detailed description of the computer program, and the input/output formats are 
presented. The results of sample calculations performed with the computer program 
are compared with experimental data. 
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2.0 -INTRODUCTION 


' A continuing problem in the development' of air breathing propulsion systems is 
the design of efficient subsonic diffusers Not only is the (engineer faced with 
building an efficient diffuser, but frequently he must tailor the exit flow within 
certain physical constraints imposed by the propulsion engine. This task lias proven 
extremely difficult in the past since the engineer has been forced to rely on 
empirical design methods based on correlations of limited experimental data. Diffuser 
performance maps based on empirical correlations have been published by Reneau 
(Ref. 1) for incompressible two-dimensional flow and by Sovran (Ref. 2) for incom- 
pressible annular flow. These performance maps generally correlate pressure rise', 
coefficient with area ratio and length and serve as useful design tools within the 
class of geometries tested. In addition, regions of stall on the performance map 
have been defined by Fox (Ref. 3 ) for two-dimensional diffusers and by Howard (Ref. 4 ) 
for annular diffusers. Thus, reasonable bounds have been set for the design of 
unstalled diffusers. Sovran (Ref. 2) has developed a correlation for pressure rise 
coefficient as a function of inlet blockage due to boundary layer displacement, thus 
adding a new parameter to the correlations. The effect of inlet Mach number has been 
examined by Runstadler (Ref. 5) for two-dimensional diffusers and performance maps 
published. Comparable performance maps for annular diffusers apparently have not 
been published. 

Even with the large effort that has been expended in the past in the development 
of empirical diffuser design methods, sufficiently generalized methods have not been 
forthcoming. Considering the vast number of geometric parameters, such as area 
ratio, length, wall curvature, and support struts, together with the physical flow 
parameters, such as Mach number, Reynolds number, swirl angle, and inlet flow distor- 
tion, it is evident that the development of such methods will be extremely time 
consuming and very costly. Clearly, the availability of accurate analytical design 
tools would enable the engineer to arrive at an efficient design for a specialized 
diffuser in a much shorter period of time at less cost while at the same time 
providing a better understanding of the diffuser flow processes. 

The development of analytical diffuser design methods has generally not kept 
pace with empirical studies. The conventional solutions, such as that used by 
Sovran (Ref. 2), divide the flow field into an irrotational free-stream flow and a 
boundary layer flow. These methods, which require an iteration between the potential 
flow pressure field and the boundary layer displacement thickness, frequently fail 
to converge when the boundary layers merge and cannot account for such problems as 
inlet swirl and flow distortion. Recently, however, Anderson introduced a new method 
(Refs. 6 and 7 ) for solving the diffuser flow problem-. This method involves 
numerically solving the equations of motion for the entire flow in the diffuser 
thereby enabling compatibility between the inviscid flow and boundary layer to be 
achieved without the need for matching a boundary layer solution to an inviscid flow 
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solution. The method provides the accurate analytical design tool required to define 
optimum diffuser designs for a wide range of applications in a short period of time 
at a minimum of cost. Good agreement in the prediction of pressure rise coefficient 
and stall has been demonstrated for both straight -wall and curved -wall unseparated 
diffusers. In addition, Anderson has demonstrated the applicability of the method 
to the treatment of the effects of inlet distortion and inlet swirl. 

This user's manual describes an extension of the method of Ref. 6, to enable the 
calculation of swirling, turbulent compressible flows in axisymmetric ducts. A 
generalization of the method to include treatment of arbitrary distributions of heat 
and mass transfer at the boundaries and the effects of struts or blades in the 
diffuser passage is also described. ; 
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3 .0 " analysis 


The present analysis represents an extension to compressible flow of the analysis 
presented in Refs. 6 and 7 for computing incompressible swirling flows in axisymmetric 
ducts. The unique feature of the analysis is that a conceptual division of the flow 
into an inviscid, irrotational flow and boundary layer is not made but rather the 
analysis solves for the entire flow across the duct at each streamwise station. The 
problems of matching solutions for the inviscid flow and boundary layer are thereby 
eliminated and no difficulties exist when the boundary layers on the duct walls merge . 
The analysis supposes that the streamlines of the actual flow through the duct will 
not be appreciably different from the potential flow streamlines; an assumption which 
from simple physical considerations must remain valid for internal duct flows with 
only small regions of separation. Thus it is possible to make boundary layer -type 
approximations, such as supposing the transverse velocity is small, in a coordinate 
system based upon the stream function and velocity potential of the potential flow. 

In this manner, the viscous effects can be treated as a perturbation upon the inviscid 
field and, as such, governed by a parabolic partial differential equation. The 
elliptic properties of the flow (i.e., downstream conditions affecting the upstream 
flow), are retained through the choice of the potential flow streamlines as the 
coordinate system for the parabolic partial differential equation. Included in the 
present analysis is the capability for treating the effect of struts, inlet guide 
vanes, and flow straightening vanes, as well as the effect of distributed wall bleed. 

A description of the various components of the analysis is presented in the 
following subsections. 
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3.1 Streamline Coordinate System 


The orthogonal streamline coordinate system used in the present analysis is a 
simplified version of the coordinate system described in Ref. 7. The coordinate 
system of Ref. 7 was obtained from the solution of the plane potential flow through 
the duct in question where the stream function formed the normal coordinate and 
the velocity potential S formed the streamwise coordinate. Rotation of this solution 
about the axis of symmetry provides the third coordinate (see Fig. 3.1.1')* As 
demonstrated in Ref. 7, use of this coordinate system provides a means of reducing 
the general equations of motion for the diffuser from an elliptic set of equations 
to a parabolic set of equations by making the well-justified assumption, for unstalled 
diffusers, that the flow normal to the streamwise coordinate is small. A second 
advantage is that since the equations of motion are written in generalized curvi- 
linear coordinates, more general curved wall duct geometries can be treated. 


The stream function n and velocity potential S used to generate this potential 
flow coordinate system satisfy Laplace's equation, that is 


i^ + i5u o 


dFT 


dl c 


(3.1.1) 




(3.1.2) 


In addition to the wall boundary condition, n =0 along the inner wall and n - I 
along the outer wall, the condition of no normal pressure gradients at the inlet and 
exit of the duct [dv/du = Q ) must be satisfied where V is the potential flow velocity 
or metric scale coefficient. Thus from orthogonality conditions, the angle that the 
potential flow streamline makes with the axis of symmetry (see Fig. 3.1.1) becomes 

COS 6= vg-v-ff (3.1.3) 

an os 


sinQ = v 


<3R 

as 



(3.1.4) 
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The curvature of the coordinates are. given by 


_L - _ ^y. 

P s ' d% 


(3.1.5) 




av 

as 


( 3 . 1 . 6 ) 


Distances along the three coordinates are determined from the metric scale coefficients. 




( 3 . 1 . 7 ) 


x =/dJ. (3.1.8) 

R <f> = f Rd$ (3.1.9) 


It should be noted that this coordinate system is not the only orthogonal 
coordinate system that can be constructed from the plane potential flow solution. 

The solution is made unique only with the boundary condition dV/dn = 0 at the inlet 
and exit of the duct. This boundary condition, as shown in Ref. 7> implies a zero 
normal static pressure gradient at the inlet and exit of the duct in the absence of 
swirl. Ducts with normal pressure gradients at the inlet and exit can be treated 
by extending the duct with the required curvature. 

Because of the requirement that the n and S coordinates be accurately determined 
up to third derivatives, extreme care must be exercised in obtaining the potential 
flow solution by numerically solving Eqs. (3.1,1) and (3.1.2) as described in Ref. 7. 
In addition, such a numerical solution requires significant computer time. 
Accordingly, for the present analysis a method of approximating the potential flow 
coordinate system which yields accurate coordinate derivatives and saves computer 
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time was used. This approximate method' has proven adequate providing the .diffuser- . 
wall does not have rapid changes in wall curvature or large diffusion angles. Since 
these conditions would generally lead to separation, which the analysis is presently 
not capable of handling, use of this approximate method represents little additional 
sacrifice in generality. The approximate potential flow coordinate system is obtained 
as follows. 

Assuming that the curvatures are small, it is possible to construct a circular 
arc of radius (see Fig. 3*1*2), such that it satisfies the orthogonality relations, 
Eqs. (3*1*3) and (3*1*4), on the duct walls and mean line. Then from Eq. (3*1*5) 

~fa~ = (V T "V H ) = ' (3*1*10) 


and, ' . ■_ 


v = + <v T -v H )<n-,/2) = v„--^- 

*3 


(3*1*11) 


The radius of the circular arc, from Eq. (3*1*6), is given by 


A. .. JVfi i 

- ds 


(3*1*12) 


where second order terms in l/yo^ are neglected-. The distances along the arc, which 
approximates ;a potential line is given by 






(3*1*13) 


and the distance along the chord in terms of the diffusion angle $ is, 

' - - i i • ■ : . ' 

Y = 2 J sin #/2 .a Y - 0 (0 3 ) ' • (3*1.14) 

r s 

• * . C * : 
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Hence for small diffusion angles, Q « I , the distance along the arc can be 
represented by the distance along the chord line. The potential line must now 
satisfy the orthogonality relations. On the wall, from Eqs. (-3.1.3) and (3.1.8), 

J_ _ 1 dz dx 

V ' a/, + (dR ) 2 dS = ds (3.1.15) 

dt 


Hence, for two adjacent potential lines 


V H ax t 
V T = ax h 


(3.1.16) 


and from Eq. (3.1.7) 


Yfci Xy z I 


(3.1.17) 


The metric scale coefficient can be found from Eqs. (3.1.16) and (3.1.17), once the 
potential lines are determined. These potential lines are described by chord lines 
(see Fig; 3.1.2) which are perpendicular to the mean line of the duct 

Rm = * R h + R t>/2 ! (3.1.18) 


satisfying the relation 

dR M = J_ r^H dR T ] = . h i 1 H (3.1.19) 

dz 2 l dz dZ J - r t - Rh 


In summary, the streamline coordinate system is constructed by first locating 
the chord lines according to Eq. (3.1.19)> calculating the metric scale coefficients 
along the wall from Eqs. (3.1.16), (3.1.17), and, finally, determining the curvatures 
from Eqs. (3-1.10) and (3.1.12). 
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FIG. 3.1.1rORTHOGONAL STREAMLINE COORDINATE SYSTEM. 
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FIG. 3.1.2-CONSTRUCTION OF STREAMLINE COORDINATES. 
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3.2 Equations of Motion 


The equations of motion for a compressible flow in a general orthogonal curvi- 
linear coordinate system have been derived by Pai'(Ref. 8). Substitution of the 
metric scale coefficients ( I / V , 1 / V, R ) into the equations derived in Ref. 8 produce 
the following equations of motion in the streamline coordinate system described in 
Section 3.1* r 


Continuity Equation 


dp 

dt 


. V 2 f d i r P u s \ V 
+ ~ l dS 1 V 1 + 


d i r P u n 
dn V 




<5 cfi 


0 (3.2.1) 


S Momentum Equation 
d u, 


du c 


ST + ^ u s V 1# + ^ u n V dn 


^Hs + pu$ du* p u n 
* ■# ■ ~ 


_ £^± sin 9 + Wj" 


VdP v‘ 

as + r 


dn 


( rr sn. ) + — ( — A ) 1 _ 

[ V ' d<b v 2 J 


a 

$4 ) ' v 


■ nn 


[~k< 

T 4>4> 


rr 


ss 


sin 6 + 


■ sn 


( 3 . 2 . 2 ) 


n Momentum Equation 


du 


at 


— + ^u s V 


du n 

+ v __ + r 


du n pu$ du n pu s ‘ 


d4> 


P°<t> u _ u n 

f- cos 8 + P — = - v 

* H n 

d i r T nn > d , T nd> , l 

* ~dn <— > + -if 1 V >J 


dP v _ r 

In + -r [ 


jL_ (Jli 

as v 


£LS ) 


(3.2.3) 


L ss 


r 4>4> 


T cos 9 + 


sn 
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4> Momentum Equation 


+ ^ 


cos e+ sine = -*-£ + ^[£<^> 


+ ^r<^» + ^ ( ^>] + ^- 8 + I ^ si " e 

Energy Equation 

_ f ai , ,. .. dl , „ v . dl +Ji± dl \ 

pT { + Us v __ + UnV +-^1 j 


li r <3 , rc ^n , , a , r Qn > , a ,_SA. . l ;_ _ 

r l Is ( “7“ } a^ ( “v“ } + a^ ( v 2 'J + 2 i T ™ *nn 


+ r ss * S 3 + T <p<f> & <p<p} + T ns *ns+ *<£n + t</>s *<f>s 


Where the stress, strain, and heat flux vectors are given by- 
Strain 


au , 


av 


"2 1 = v as " Un an 


& nn . au n ,, a V 

2 = v an s as 


e J± . x au <£ . „ V ar v ar 

2 - t “a^r + u s— is- + Un ^- aT 


a , u< £ . i o , , 

*n<£ = r V ( — ) + J- y -gy (u n V) 


(3.2.4) 

(3.2.5) 

( 3 . 2 . 6 ) 

(3.2.7) 

( 3 . 2 . 8 ) 

(3.2.9) 
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{ ^ v,+ rv -k i! r- ) 

«»■ aF (u nV)+-|-(u s v) (3.2.11) 

Stress 


r ss = H- *ss - P Os o s ' - -§-M v 

(3.2.12) 

T nn = H-enn-P u nOn - -y* 

(3.2.13) 

— . — j- 2 

*00 = H- e 4><f> -yoo'^u^ y AtV 

(3.2.14) 

*h0 = h- 8 n<p ~ P u n u £. 

(3-2.15) 

^s0 = H- g s4> ’ /> u s u 0 

( 3 . 2 . 16 ) 

T sn = /^ *sn ‘ /» 0 s u'„ 

(3.2.17) 


( 3 . 2 . 10 ) 


where the stresses include the Reynolds stresses and the dilitation V is given by 




(3.2.18) 
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Heat Flux 


H-E d 

<in* TfeT v ^r(Cp T > 
V t£ e v i? (c p t) 


(3.2.19) 

( 3 . 2 . 20 ) 
( 3 . 2 . 21 ) 


Since flow losses are of special interest in diffusers, the energy equation, in the 
present analysis, was written in terms of the entropy, that is 




DI r DT - -DP - 
Df " P Cp Dt Dt 


Q + i> 


( 3 . 2 . 22 ) 


It can be seen from Eq. (3.2.22) that the change in entropy is the small difference 
between two large terms and involves only dissipation and heat transfer. The 
thermodynamic variables can be recovered from the entropy by the relation 


I f = Cp In (-y-) - In (-^“) 


(3.2.23) 


With struts in steady flow, circumferential variations of the flow variables 
exist. In the present analysis, these variations were averaged according to the 
general relation 


f(n.s) 





f (n,s,c£)d<£ 



lU 


(3.2.24) 



Writing the dependent variables in terms of the mean f and deviation f as (see' ; 
Fig. 3.2.1) 


f(n,s,<p) = f (n,s) + f (n,s,<£) 


(3.2.25) 


and then integrating 


r* 2 li 
\ d - a 


d< £ ' da t R ^ ] + ^ 


* Jl 

da 


F d< t>2 
2 5 a 


(3.2.26) 


L 


$2 df 


d< £ = f 2 f I 


where for the case of no struts 


f 2 = f\ 


(3-2.27) 


(3.2.28) 


These above rules of averaging are similar to those used in obtain' ag the 
Reynolds stresses from the Navier-Stokes equations (Ref. 9) and are appi.Led to the 
present equations of motion accordingly. With the assumption that 


f'/f «l 


(3.2.29) 


all double and triple correlations can be neglected except those terms involving the 
stresses in the strut boundary layers and the stresses produced by the mean flow which 
produce the annulus boundary layer. Since the equations of motion are written in a 
coordinate system approximating the real streamlines, the additional simplifying 
assumption can be made that 
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Up /Ug <<: I 


(3.2.30) 


and that only the shear stress terms Z ns and and strain terms E ns and E n< ^need 

be retained . 

Application of the above averaging procedures and simplifications results in 
the following basic equations of motion in which H s , Hj are terms identified as 
strut forces and <$> s as the dissipation function for- the strut boundary layer. 

Continuity Equation 


J_ (d± 
d7l Us 


Ts - v p 


v as 


= pu s 


(3.2.31) 


S Momentum Equation 


p au s v au s a* 1 an 
s as ~ g aw as 7 K 5 " ds 


jL [ G ^ ns ] _ ^ns av + H s u* a r 
g aw L v J v a» v *r a s 


r as 


(3.2.32) 


Momentum Equation 


_ v ay _ v a 

as g ah as "S" aw - 


dn as 


v _a_ 

IT dll 


(3.2.33) 


dR H <£ p U s U ^> a R 

r aw v “ r as 
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71 Momentum Equation 


an 


y m 2 a% \ as 

T 1 av ..2 . i aR .. zH 
L v d% Us R dV. U ^J 


Energy Equation 


'+P 




ap 

as 


a 2 R 


LR a % as 


av 

au s 

— 9 

PU<£ 

3R dU 4> 

dn 

a s 

c. 

R 

d% as 

1 

*d 

c 

m ro 
1 1 

\ a 2 v 
.v 372 as 

1 av 
v 2 d % 


1 

aR 

aR 



R 2 

a 7i 

as _ 


1 

If- 

-7 

V 2 

a r go 1 

v® 1 a 

- 1 

G dTL L V J 

In 

cfi ^ n 

'</> + 

V 

]} 


( 3 . 2 . 3 *+) 


( 3 . 2 . 35 ) 


where the stream function is expressed as 


av// _ _g_ 
an V 


PU S 


( 3 . 2 . 36 ) 


d'p 

as 



( 3 . 2 . 37 ) 


and the gap G is expressed as 


G = R(<£ 2 -<£,) 


( 3 . 2 . 38 ) 


The continuity equation and the s and </> momentum equation and energy equation follow 
directly from the equations of motion in the streamline coordinate system. Equation 
(3.2.3*+) was obtained by differentiation of the n momentum equation given by 
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(3.2.39) 


an 

3 % 


r 


2 _E. iy. 2 p. _a_R 

r V v Us r dt 


u 


%} 


Equation (3.2.39) illustrates the significance of the approximations used in the 
streamline coordinate system; namely that the pressure field and particularly the 
normal' pressure gradient is determined by the radius of curvature of the streamline 
coordinate system rather than the radius of curvature of the real streamlines. 

To complete the set of equations for the number of unknowns the equation of state, 
definition of entropy, and the definition of the stress strain relations and heat 
flux given below, are used. (The equation of state and entropy equation, (i.e., Eq. 
(3.2.23)), are differentiated in the s direction for direct application of the explicit 
numerical integration method.) 

Equation of State 


an _ p a_© 
as r a s 



= o 


(3-2.40) 


Entropy Equation 


ai x i a® i an 

as a-i ©as n as 


Stress, Strain Components, Heat Flux 


_ a u s a v 
Ens= v In + Us Tr 


F - V _ w U< £ d R 

E n<£ v a 5? v R 3 91 


'ns 


H-r 


E 


ns 


(3.2.41) 


(3.2.42) 


(3.2.43) 


(3.2.44) 
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(3.2.45) 




E n^> : 
Mr 


Q 


I 

n r p re 




(3.2.46) 


In addition the following thermodynamic relations such as the isentropic flow 
relations, and definitions of engineering parameters are used. 

Thermodynamic Relations 

n = P© (3.2.47) 


i = |n ®~ ln n 

M - M r U/C 


c - y© 


np 

n 


[• 


+ 


7-1 

2 


M 2 ]/ -1 


©o 
: © 



M 


Z 


(3-2.48) 

(3.2.49) 

(3.2.50) 

(3.2.51) 

(3.2.52) 
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Additional Relations 


= 2Z 
* Poo Um 


(wall friction coefficient) (3.2.53) 


Rc = «^/©cd l a M co 


(recovery factor) 


(3.2.54) 


Q 


w 


P CO U 0O (@ CO' @ A 


(Stanton number) 


(3.2.55) 


S t = Cj/( 2S a ) 


(Reynolds analogy factor) (3.2.56) 


These equations of motion form a complete. set of equations for the solution of a 
flow problem provided a priori knowledge is obtained for the blade forces H g , , 

and for the blade losses . This information is obtained from the , blade section 
lift and drag coefficients and is described in Section 3*3. 
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3.3 Blade Element Forces and Losses 

The "blade element forces are obtained by collecting those terms from the 
equations of motion derived in Section 3.2 that do not cancel out at the strut 
surface. These terms are 


r V 


d4 > ; 


h z — 

+ ( T sn 2 dn 


2 “as" 

d<t> 


2. . 


fn 


rv r d4> 
— [ <p - — 


( T r 


2 dn 
te 2 


nn2 ,j n 


rv 
d<f> 2 
r n<*> 2 “aT” 


a<£ 



d 4 > 2 


d<t>! 

P| as 

(T SS2 

as 

r SSI 

as 


d<t> | 


T te 2 - T 

tei ■ 

i 

T sn i 

an 

’ ) + - 

r V 


1 

64 >i 

) - 


CVJ 


te\ 

p i an 

( r sri 2 

as 

r sn i 

as 


d<t> 


CVJ 

-e 

e 



T nn i 

an 

- ) + 

rV 

]* 

P, , 


d<t> 2 


it. 


— -(T 

S<£2 

as 

T S<£l 

as 

) 


te\ ^ 

T te I ~drT 5 + Fv J 


(3.3.1) 


(3.3.2) 


(3.3.3) 


It can be seen that these terms contain all the skin friction forces and the normal 
forces produced by the static pressure at the strut surface. 

To demonstrate that the above terms are precisely the forces acting per unit 
volume on the flow, consider an element of volume (see Fig. 3-3.1) given by, 


dS dn 
dV = g — — 


(3 .3 -*0 


and an element of area on surface (l) by 


dA 


dS dn 

~V ~V 


V I + (rV V ! + (rV -“TT ") 2 

dn dS 


(3.3.5) 


such that the direction cosines of the blade surface vector, "/p normal to the blade 
surface (1) are 
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V is 


X '"' '"• = J 


d<j>> 

- rV -^ 

v / i + (rv^J-) 2 + (rV 

dn 
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80, 
~ rv an 



-v d c6 1 z %/ d4> i 72 

l+(rv— -i-) l + (ry-^i-) 


dn 


dS 


Then the force acting on an area dA, for surface (l) is given by 
^ * 1 

/, = {[(- P + r ss ) i s + r sn i n + r s £ i^j • X, } i £ 

+ f [ T ns 's + ( " P + -' Tnn) in + *"n<£ i$] * X, } in 

+ {^s >s + i n + i<£ J • y, } i<£ 


(3.3.6) 


(3.3.7) 


(3.3.8) 


(3.3.9) 


This produces a force per unit volume acting on the flow from both surfaces 

— ^2 dA 2 - /j dA, __ _ — 

f ~ (jY f s 1 s ¥ / n 1 n 1 <jf> 


( 3 . 3 . 10 ) 


■where / s ,f n , / cf> are precisely the terms defined by Eqs. (3-3.1) through (3-3.3). 
The total blade force is obtained by integrating over the volume to obtain 

... . .. 7'/7a*S- nf- (3.3.H) 

¥ 
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(3.3.12) 


Hence the blade force per unit span is given by 



and the blade force per volume is related to the section blade force by 


f A1 

1 ' 9 dS 


(3.3.13) 


In the present analysis provision is made for obtaining the blade force from 
either the section lift or drag coefficients of isolated airfoils (the case for 
struts with low solidity) or from loss coefficient and turning angle correlations 
(e.g., Refs. 10, 11, and 12) (the case of closely spaced struts or blades). Before 
the blade force can be obtained by either method it is first necessary to determine 
the streamwise static pressure differentiation (back pressure) across the struts, 
which can be significant, particularly if the flow is swirling. Since the back 
pressure depends upon knowledge of the blade forces, (however, an iteration between 
back pressure and blade forces) is required. This iterative solution is obtained in 
the following manner using methods developed in Ref. 12. 

To develop the required .iterative solution it is first necessary to develop 
relations for determining the blade forces f $ and f<£ and the blade dissipation <J> B 
(or, in dimensionless form, H s > and <J> B J in terms of the flow variables 

entering (station 2) and leaving (station 3) the blade row (see Fig. 3.3.2). 
Considering the flow in a blade element strip Ay =An/V , as shown in Fig. 3.3.2, 
and assuming that U n is small, and the walls are adiabatic, the S and momentum 
equation, with the additional assumption that the shear stresses can be neglected, 
and the energy equation in dimensionless form (Eqs. (3.2.32), (3.2.33), and (3*2.35)) 
become 


H; 


v an 

y M j? 3 S 


+ pu s v 


du s 

as 


(3.3.1*+) 


H<£ = PU S V 


du-0 

as 


(3.3.15) 
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(3.3.16) 


©PUg <3 (A I q) 

= r v a q 

y M r d S 


Allowing for small spanwise flows so that the back pressure can be in equilibrium, 
U s3 need not equal Us 2 > and the solution for the blade element force will be a 
function of the mass flow ratio. Making the approximation 


, (PU S ) 2 + (PU S ) 3 

(PU S ) A - „ 


(3-3.17) 


the forces H s and H<£ may be .integrated over the chord using Eq. (3.3.13)> to yield 

u,j-u«]} : (3.3-19) 


H,’. S./H, f 


y M r 


= G B 


/ 


(PU.) , + (PU a ), 


[ U (£3 _U c£2 (3*3.19) 


Differentiating Eqs. (3.3*18) and (3-3.19) yields 'for the local blade force 



(3.3.20) 


H<£ 



(3.3.21) 


The blade dissipation function $ s can be determined from the entropy change 
through the blade row. Defining the loss coefficient as 


Z 


^02 03 

B n 02 ■ n 2 


(3.3.22) 
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the entropy rise becomes, using the isentropic flow relations, Eqs. (3.2,51) and 
(3.2.52), as well as ' the' assumption that the total temperature is conserved since 
no work is done, 

(3.3.23) 

From Eqs. (3.3.18), (3.3.19) and (3. 3.23)H S , H<£ , and $ B or 3 S , 3 , A I B , are 
known in terms of the flow variables upstream and downstream of the strut. Having 
developed relations between the blade forces and dissipation function in terms of the 
upstream and downstream flow variables, it is now necessary to develop additional 
relations for computing the downstream flow variables given the upstream flow. 

Assuming that the blade force in the annulus boundary layer is the same as the 
blade force at a distance from the wall equal to the displacement thickness of the 
boundary layer, the flow variables at station (2) and station (3) can be calculated 
as though the flow were inviscid, This is accomplished in the present analysis by 
setting all the flow variables in the annulus boundary layer between the wall and 
the boundary layer displacement height equal to their value at the displacement 
height, that is 

'■ (Y) '= c/>(8*) ' for 0 < Y < 8 * h : . 

$ (Y) = 4 > (Y) 

$ (Y) - (Y-8*) 

The solution of this inviscid flow problem requires that the flow satisfy the . 
equations of motion (Eqs. (3.2.31), (3.2.32), (3.2.33), (3*2.35), and (3.2.39)) with 
the viscous terms deleted.' ,The S and <f> momentum equations (Eqs. (3.2.32) and 
(3.2.33)) are satisfied by introducing lift and drag coefficient correlations, or 
turning angle and loss coefficient correlations depending on which is known. 

If C L and c 0 are known then and are determined by a simple resolution 
of forces 


(3.2.25) 


^ 4 >'~ G B^ 2 [J z a [ -c L cos ° 2 +C 0 $in a 2 ] (3.2.26) 


A A 2 

3rr = - G„P, U- <T 


[-C L sino 2 + C D co$Q 2 


8h< y < Y t -8 h 
Y T -8y < r < y t 
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and Eq. (3-3.18) and Eq. (3-3-19) may be used to solve for u s3 andu^, 3 - The energy- 
equation (Eq. (3-2-35)) in this case reduces simply to 


@ 03 -' @02 


(3-3.27)- 


since no work is done by the blades. The remaining equations to be solved are the 
continuity equation, Eq. (3-2.36), and the normal pressure gradient equation, 

Eq; (3-2.39) > which are expressed as 


_ GPu s (3.3.28) 
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(3-3-29) 


If the exit flow angle a and loss coefficient Zg are known, Eqs . (3.3.28) and 
( 3 - 3 - 29 ) may be reduced tp 



(3-3-30) 
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2H 
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]} 


(3-3-31) 


where Eq. (3.3-22) is used to determine Il 0 from Zg. 

It is evident from the above discussion that for both the case where C L and C D 
are specified and the case where a and Zg are specified, there are two differential 
equations relating the blade foces and dissipation to the back pressure; that is, 

Eqs. (3-3-28) and (3-3-29), or Eqs. (3-3-30) and (3-3-31). These equations are solved 
iteratively using the boundary conditions. 

(3-3-32) 
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3.4 Turbulence Model 


The 'turbulence model used in the present analysis is the same as that used in 
Ref. i'3 with the inner layer region governed by Prandtl’s mixing length theory 
(Ref. 14) with a Van Driest damping factor and the outer layer region. governed by an 
eddy viscosity model proposed by Clauser (Ref. 15). The inner. layer solution is 
obtained analytically, using the method of Spalding (Ref. l6)'. and the_outer layer 
solution is obtained numerically. • ; .. .. 

Outer Laver Solution . . ■ ' . , 

For ^Clauser 's eddy viscosity model, we have from Ref. 15 

/*E = X/° u co S * .... . .. (3.4.1) 


where Uco i- s the velocity at the edge of the boundary layer and 8* is the boundary 
layer displacement thickness. In dimensionless form the eddy viscosity can be written 


= xNrPucoA*. 


(3.4.2) 


In the present analysis Uco 4 s defined as the maximum velocity in the duct and the 
average velocity U and /Of are defined as 

2, UoD A*= /(Uco-u) ^ =(U-U)Y T " J (3.4.3) 

0 


Thus the outer layer eddy viscosity is (see Ref. 7), 

7 f = XN R P ao| I (^-0) (3.4.4) 


where for thick boundary layers or for flows with radial total pressure distortion 
the mixing length is of the order of the duct width. For thin boundary layers, 

Eq. (3.4.4) reduces to Clauser' s law of the wake model. 
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Inner Layer Solution 


In the region near the wall, the independent variables and their derivatives 
change rapidly and numerical solution is difficult. except with a large number of . 
closely spaced mesh points. Using the numerical method chosen for the present 
analysis, the step size in the streamwise direction must be very small to obtain a 
stable numerical solution with a large number of closely, spaced mesh points; thereby 
resulting in large computational time. However, Spalding (Ref. 16) has shown that 
the solution for the flow near the wall can be computed from analytical relations 
with reasonable accuracy. -Therefore, Spalding's analytical method was used for the ■ 
inner layer calculation, to enable larger step sizes to be employed in computing the 
diffuser flow. This method, however, may not work well for very high Mach number 
flows, or when the temperature and momentum displacement thicknesses are very 
different. 

According to Ref. 16 the equations of motion for the inner layer in universal 
turbulent coordinates are given by, 


Me ^ du_ _ r » = j + p+ y + + m* u + • . - 


Mw / dY 


... (3.4.5) 


'_£!-) JL«; = 0* » , + M + H + 
Pre Mw / dY 


(3.4.6) 


i Me \ d 

~ ” '“A 


Pre Mw dY 




( 3 . 4 . 7 ) 
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(3-4.8) 
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^ (3.4.9) 
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and for the stream function 



2«* = u* 

dY* 

(3-4.10) 

For large values of Y + , Eqs. (3-4.5) through (3-4.9) simplify to 



du + . Vt" + 
dY*- KY + 

(3-4.il) 

■ ■ - 

dH + . p rt q+ 

du* • t * 

(3-4.12) 


dH A - (p )u+ Pr T M+H A 

dU + P ^T ,)U + T * 

1 .'•••• ..-i 

(3-4.13) 

which yield 

‘ f \ 

the following analytical solutions for the cases of interest. 


CASE I m+ = 

P = 0 (Flat -Plate, No Bleed) - 



U* = y In (C 5 Y + ) 

..(3-4.14) 


w+ = Y [r ln(C 5 Y + )-Y + ] +C 6 

(3-4.15) 


u +2 - 

• = (P RT - 0 T +X 2 

(3-4.16) 

6 

H + = U + + C 3 

(3-4.17) 
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CASE II P*~ 0 » 0 (Flat -Plate, With Bleed) 


u* = -£-ln(c,Y*) + [ln(c,Y^ 


(3.4.18) 


W + = ^[ln(C 1 Y + )-l] + ^- 2 {v + [ln(C l Y + )] 2 -2Y + [ln(C l Y + )-l]} + C 4 (M + ) (3.4.19) 


Ha = 


( P RT 


U + 2 (|+mV) P rT -[| + P rt M + U + + P rt (F{ ?t - l)(M + U + ) Z /2 


('-PrtKS-Pr^M 


Uc 2 (I+M 


V) Prt 


(3.4.20) 


H + = -^{i-exp [in (i+m\0 + C 3 M* ]} 


(3.4.21) 


CASE III M* r 0 , P* > 0 (Pressure Gradient, No Bleed) 


U = ^ In (C 5 


mi- [yr^ -q-±*[i±* 2±g±rz-] (3.4.22) 


w + = [y + In (c 5 Y + ) -vj + 1 [-^(i + P + Y + ) 3/2 -Y + - -~r ] 


- ~ { Y + ln 


2 / l+/T + P^ Y* ■ Vh-p^y*- 


(3.4.23) 


Y — I 1 
3+ J 


+ c. 
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(3A.24) 


H A + = (P RT -0^. + C 2 

i < 

H* = U* + C 3 . , (3.4.25) 


Note that for Cases II and III, in the limits as M + -*-0 and P + ^-0, Case I is obtained. 
The constants of integration in Eqs. (3-4.14) through (3.4.25) are obtained by 
numerically integrating Eqs. (3.4.5) through (3.4.7) out to a value of Y + and eval- 
uating the constants by matching with the analytical solutions. It should be noted 
that since analytical solutions cannot be found for the general case, Case III is 
used for flows with pressure gradient and bleed. i ; 

The analytical solution is used to compute the flow in the inner layer which is 
matched to the numerical solution for the outer flow at some streamline near the 
wall. This matching is accomplished by iteratively finding the friction velocity U* 
such that Eq. (3-4.18) or Eq. (3.4.22) is satisfied on the streamline where the flow 
is matched. Once the friction velocity is found the adiabatic wall temperature may 
be determined by matching the temperatures at the matching streamline using 
Eq. (3.4.20) or Eq. (3.4.24). If the wall temperature is specified, Eq. (3-4.25) is 
used to determine the heat flux. 
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3*5 Classification of Equations 


In fluid mechanics, the equations of motion are generally classified as elliptic 
for inviscid subsonic flow, hyperbolic for inviscid supersonic flow, and parabolic 
for boundary layer flows. The classification of the equations being solved in the 
present analysis is required to establish the proper boundary conditions and to 
evaluate the general properties of the solutions. 

The equations of motion to be examined are Eqs. (3*2.31) > (3*2.32), (3*2*33) > 
.(3*2.35), and (3*2.39) with the stresses and heat fluxes substituted appropriately. 
Then we have five equations , continuity, three momentum equations, and the energy 
equation for the five variables, U s , U<£, U n , II , © * If a is the characteristic . 
surface, then according. to Ref. 17, the following determinant vanishes on the 
characteristic surface. * 
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• . .0 
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C 3 a, n 

C 3 a, n 

0 


C 4 a 2 


where the subscript notation is used to denote partial derivatives and 


* = U s £2,5 + U s a, n 


(3*5*2) 



C 2 = l/( r M r 2 ) 


c 5 = yMr 2 V 2 ^/P 

I 7 y (Pl\ 

4 " It y-\ Nr p R \ Mr / 


(3*5*3) 
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Evaluating the determinate, we have 



C 4 cfPa,® =0 


(3.5.4) 


Hence all roots are real and equal and the equations of motion derived here are 
parabolic . 

Since the equations of motion are parabolic, the proper boundary conditions can 
be prescribed. Along the inner and outer walls 


u sw = w = 0 
u nw = M w / Pw 


(3.5.5) 


For the energy equation we have either adiabatic walls or prescribed temperature 


Q w = 0 or-©=© w '•••. (3.5.6) 


Finally, the flow entering the duct U s (n), U^tn), II ( n ) , ©( n ) must be specified. 

The exit flow conditions cannot be specified but must be determined by the solution . 


The subcharacteristics (Ref. 18) are determined by the equations of motion 
minus the diffusion. terms and the resulting determinate is given by, 


D 


a.s 


0 


a, n i///n 


Pi// 0 0 

0 Pi// 0 

• 0 . . 0 _ [p 1 //.] 

000 


62 G) s 0 

0 v : 0 

■ • 1 ■■ 

C 2 C2n ' - - 0 

- ^ yrj- 


(3.5.7) 
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where the boundary layer assumption has not been made and the term in the bracket is 
the term that will vanish when the boundary layer assumption is made (i.e., it contains 
the acceleration terms in the momentum equation). Evaluating the determinate 


D = 


\j/3 y 

y M? ~y 


7 j- a,s + 5 2 M s M n a lS a,n+ (EM^-Oa^} (3-5.8) 


where 


e =• [p^/] / (P\p) 


(3.5.9) 


When the boundary layer assumption is made ( £ = 0), the characteristics are real 
and equal. If £ = 1, the characteristics are imaginary for M <1 and real and 
distinct for M > 1. Thus, it is seen that if the boundary layer assumption is not 
made, the subcharacteristics have the properties that are usually associated with 
fluid dynamic flow regimes. It is only the boundary layer assumption which reduces 
the equations to a parabolic system. 

Although the equations of motion are parabolic for subsonic flows, the solutions 
can be expected to have some properties of elliptic flows since the coordinates were 
obtained from the solution to an elliptic differential equation (Laplace's equation). 
Thus, normal pressure gradients required to turn the flow are a function of the 
curvature of the streamline coordinates and the calculated flow starts to turn 
upstream of a corner. However,- the upstream flow cannot start to turn due to 
downstream flow separation since the streamline coordinates do not have this 
information. 

For supersonic flow, the equations of motion are still parabolic, but the 
inviscid portion of the flow should have properties of solutions of hyperbolic 
equations. Discontinuities, such as shock waves may exist, and turning of the flow 
at a corner due to. an expansion fan may exist in hyperbolic flows. These flows 
cannot be described well by the present analysis because the streamline coordinate 
system currently incorporated into the analysis does not accurately describe the 
curvature of the real streamlines. However, in duct flows for which the turning is 
more gradual, and for which the Mach number is close to 1.0, reasonable solutions 
may be possible with the present analysis since the Mach lines are nearly normal to 
the flow and may be approximated by the potential lines, which are the characteristics 
of the parabolic equation. In principle, it may be possible to solve the problem 
for supersonic flow with expansion waves and weak shocks provided an orthogonal 
coordinate system can be constructed from a supersonic inviscid solution. It should 
be noted, however, that for strong shock-boundary-layer interaction, the method would 

require several iterations because the shock location would depend on the boundary 
layer growth. 
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3.6 Method of Solution 


The equations of motion, Eqs. (3.2.31) through (3.2.35) and Eqs. (3.2. 40) and 
(3.2.4l) are solved by the explicit numerical method of Ref. 19 . This explicit 
method was chosen on the basis of its success on the analogous incompressible flow 
problem described in Ref. 7 after other numerical methods had been tried. With this 
method, the solution for the dependent variables and their normal derivatives are 
presumed known from previous calculation or input data up to station S. After some 
manipulation, these seven equations may be reduced to two equations of the form 


JL [ ] - 

d% \ as 




an 

as 
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(3.6.1; 


±( m - f 

39l \ as / n Las • as 
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(3.6.2) 
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These -equations 'are then linear in the unknown dependent variables <3 \J//d:Sand -dlT/fdS ; •: 
Thus- they can -be -'thought of- as ; two linear 'ordinary differential. 'equations across ethe- 
duct which -may " l be 'solved with the boundary conditions iy -o y .ra: . •• , 
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(3.6.3) 


= ii 


as 7n = i 
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(3.6.4) 


Since Eqs. (3.6.1) and (3.6.2) are linear, they may be solved numerically. by- -g-.A: 7. > 
calculating the homogeneous and particular solutions with the boundary conditions 


O 



(3.6.5) 




'dUT' 

K d s ; / h 




and then adding these solutions to obtain the general solution, that is 


38 


( 3 . 6 . 6 ) 


. fJJ A / d_±\ , 

as \ds /p \8S/h 



The constant of integration C, is obtained by satisfying Eq. (3.6.5) at the tip wall, 
that is 


Mr . (d±\ . (±±\ / dJA 

V T \ds h* I V<JS )p L \ds J h 


( 3 . 6 . 8 ) 


From Eqs. (3.6.5) and(3.6.7) it can be seen that the integration constant is equal 
to the streamwise; pressure gradient on the hub wall.- (it is important to note from, 
the above discussion that in the present method of solution the streamwise pressure 
gradient is obtained by satisfying continuity without the need for iteration.) 

In the numerical solution of Eqs. (3.6.1) and (3.6.2) normal derivatives are 
distorted by the transformation 


X(n) = drj/dTl 


(3.6.9) 


For any dependent variable <f>(K), the normal derivatives are obtained using central 
differences, that is / . 


» 

d<f>\ 0(K + l)-<£(K-l) 

Jn) K 2A77 X (k) 


( 3 . 6 . 10 ) 


where Atj is equally spaced. The integration of Eqs. (3.6.1) and (3.6.2) takes 
place using the trapazoid rule, that is 
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/djp\ _ /< 3 \ A77 r . F ^ (K+I ^ I 

\as J K +r V^s j K + 2 L X(K) + XlK+D-i 

(dTt\ = / an\ _Atj r FrfK) FJ](K-H) i 

\dSj KM \dSj K 2 l X(K) X(K.+ I)J 


It is noted that F^( K+ 1) and Fj^( K + 1) are not known; however, they are linear 
functions of dt/z/dSand d r}/dS. Thus 

/d\JA / <3 1! : 

^v^VK^ vv), ^ * ( 3 . 6 . 12 ; 

V 4 n (K, + B n (K) (^) K +c n (K '(ir) K 

Hence., substituting Eq. (3.6.12) into Eq. (3.6.11), d^/dS and dri/dSat K = K + I 
are obtained from the relations 
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The homogeneous solution is then obtained by setting = Apj=0. 



Once d \j//dS and <3 Il/dSare determined across the duct, the streamwise derivatives 
of the other variables are found, and integration in the streamwise direction takes 
place using the second order formula 

n(S+ dS) = n(s-ds) + 2AS^-~ j (3.6.15) 


A special case exists for the boundary condition evaluated at the axis of 
symmetry since the equations of motion are singular at the origin. For this case, 
the solution may be expanded in a series about the origin and evaluated at a control 
surface with a small radii, that is 

' ‘ / 

■ - ‘ • u s '= u s (0) + C, R 2 

IU = C 2 R (3.6.16) 

! ■■ ; , • - * ' h •. ■ 1 ■ • ■■ . " ' ’ \ ’ :J '~' ‘ 

® = 0 ( 0 ) + C 3 R 2 

where can be evaluated at the first two mesh points. 



4.0 GENERAL DESCRIPTION OF COMPUTER PROGRAM 


This section is intended to describe the general features of the computer program 
in sufficient detail so that the program can be run successfully by the average 
intended user. The first subsection describes what problems can be solved and what 
problems cannot be solved. It also describes any special care which should be used 
in exercising the various options. The second and third subsections present a 
detailed description of the input which' is required in the operation of the computer 
program and the interpretation of the printed output. Since any complicated computer 
program may fail due to inconsistencies in the input or failure of the, theory, the 
computer program is provided with self -diagnostics which notify the user of the type 
of failure. The last subsection deals with these program diagnostics as well as 
helpful hints to correct problems which may be encountered. 
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Program Language 


4.1 General .Features of , Program 


The' computer' program was written ' and developed in FORTRAN V computer language 
for use on a UNTVAc 1108 computer. A special UARL INCLUDER code was used to generate 
a new source file to facilitate’ conversion of the program for, operation on an IBM 360 
computer. The actual conversion was done by the UARL CONCUR code which provides a . 
means for converting UNIVAC 1108 CUR tape source elements into IBM 360 card image 
form. The program was then compiled ’in FORTRAN IV and executed on an IBM 360 computer. 

Types of Fluids 

The fluid may be any compressible' gas as defined by its thermodynamic properties 
R, Cp, Cu , [A, P RL , P RT . If not otherwise specified, the gas is assumed to be air. 

The reference conditions for the gas properties must be specified at standard sea 
level conditions. 

Types of Flow Situations 

Subsonic, turbulent, swirling, or nonswirling flows in axisymmetric ducts may 
be calculated including flows with radial total pressure distortion. Two-dimensional 
flows may be calculated by constructing an annular duct in which the inner to outer 
radius approaches 1.0. When struts are included in the duct, the solution is given 
for the circumferential average of the flow contained between two adjacent struts. 

Duct Geometry Options (I0PT3) 

The flow through any axisymmetric duct may be calculated provided that the flow 
is generally in the axial direction. Duct flows normal to the axis of symmetry, or 
which reverse direction cannot be calculated due to logic limitations in Subroutine 
C00R. Ducts with sharp discontinuities, such as a step, which produce separation 
also cannot be calculated. 

Provision is made in the program to either read the duct coordinates from input 
data cards (I0PT3=2), or to calculate the duct coordinates analytically (l0PT3 = *4) 
from a few input duct shape parameters. If the duct coordinates are read from input 
cards, care should be taken that the input coordinates have sufficient smoothness to 
calculate the first and second derivatives using numerical finite -difference equa- 
tions. When the second option is used (I0PT3 > 4), the user must program his own 
calculation in Subroutine GDUCT. Sample programs (I0PT3=1, 3, 4) are given in 
Subroutine GDUCT for the users reference. For ducts with no centerbody a zero radius 
must be specified. 
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An important restriction to the computer program is that the inlet and exit 
flow must have no normal pressure gradients produced by streamline curvature although 
it may have normal pressure gradients due to swirl. Many ducts 'do r hot satisfy this 
requirement; however, these ducts can still be treated if the duct is extended. For 
curved annular ducts exhausting to' atmosphere, the exit' flow may have ’’'curvature . 

This phenomena may be simulated by extending the duct to approximate : the curvature - 
of the exit flow. • / 

Inlet Flow Options (igiPTl) ■ - " : •' £•*' 

The computer program is provided with two methods to describe- the inlet flow. ' 
When I0PT1=1, the inlet flow is calculated by prescribing the ''stagnation conditions" 
(P 0 , "I"o ) , the inlet Mach number M|, the swirl angle a] , and the boundary layer " 
parameters 8* and n which are the boundary layer displacement thickne'ss and power 
law velocity profile exponent, respectively. The core flow is then calculated from- 
isentropic flow relations, and boundary layers added using power law ’velocity profile 
relations. When stagnation conditions are not specified, the calculation assumes 
sea level conditions. , .• [ 

When I0PT1=2, the inlet flow is prescribed from input data cai<ds which specify 
the stagnation pressure P 0 , static pressure P ,' swirl angle a, and' stagnation tem- 
perature T 0 , as a' function of the fractional distance across the inlet. This data 
need not be specified at equidistant points since a linear interpolation is used. 'tb v 
specify the data at the mesh points used in the calculation. If experimental data is 
not used, care should be taken that the data is self-consistent and that it satisfies 
the radial equilibrium equation. Since the initial growth' of the boundary layer is 
sensitive to the wall shear stress, data describing the boundary layers should be 
accurately specified. When this is not possible, boundary layers may be added to 
each wall by specifying 8* and n . Special care should be exercised in using I0PT1=2 
option, with or without the feature of adding in the wall boundary layers. If the 
stress distribution across the duct is not smooth and realistic, numerical instabil- 
ities might originate in the inlet flow and grow rapidly to a point where the calcu- 
lation is terminated. This may take the form of an unrealistically early separation. 

Boundary Conditions 

Either the adiabatic wall or the heat transfer case may be calculated. The 
program assumes adiabatic walls unless the wall temperature is specified. Any wall 
temperature distribution may be specified, either on input cards when the duct 
coordinates are read, or calculated when the duct coordinates are calculated. The 
case of wall bleed may also be treated in a similar manner; wall bleed flow rate is 
zero, unless otherwise specified; - - - - ■ ' ' * - 



Force Options (I0PT2) 

Subroutine F0RCE is provided with four options (see Flow Chart 7-2. 7.1) for 
calculating the strut force. For I0PT2=1, the strut force is calculated from 
experimental data in the form of P Q , P, a, T 0 across a control surface just upstream 
and downstream of the blades. This mode of calculation should be very useful in 
evaluating the influence of struts on duct flows using experimental measurements. 

For I0PT2=4, the strut forces are again calculated from experimental data, but this 
time from exit angle, a 3 , and loss coefficient, Zb, data. Experimental data for 
turning vanes are more likely to be presented in this form and this option should be 
useful in evaluating experimental data and developing empirical coefficients. For 
I0PT2=2, the blade force is calculated from empirical relations contained in 
Subroutine CASG.-: For this option, Subroutine CASC. produces the section lift and 
drag coefficients of the strut. For I0PT2=3, the blade force is also predicted from 
Subroutine CASC. , In this case, however, CASC produces the exit turning angle, a 3 , 
■and the loss coefficient, Z B . 

Blade Shape Option (ISHAFE) 

The blade section lift and drag coefficients, or exit flow angle and loss 
coefficient, are calculated in Subroutine CASC according to the blade shape index 
I SHAPE. Since hundreds of blade shapes and turning vanes exist and cannot all be 
programmed, the user must modify Subroutine CASC for his use. However, as an aid to 
the user in developing this subroutine, two options have been programmed. For 
ISHAFE=1, the struts are NACA 5 digit series uncambered airfoils. The output for 
this option is in the form of section lift and drag coefficients. For ISHAFE=2, 
empirical relations for turning vanes are programmed. In this case the output of 
Subroutine CASC is the exit flow angle and loss coefficient. 

Failure Modes 

. ■ *. C 

In the event of failure in the calculation, the program prints an error message 
called "diagnostic". These "diagnostics" are in addition to the computer diagnostics 
and are clearly labeled as such. These "diagnostics" terminate the calculation only 
when very serious. A list of these "diagnostics” appear in Section 4.4. Included 
with this list is an identifying number for the "diagnostics", the location 
(Subroutine), and the immediate cause of the failure. Where possible, suggestions 
are made to correct the calculation. 

Debug Options (IDBGN) H . 

r V 

Auxiliary printout which was originally used to debug the computer program is 
available to the user by setting the appropriate IDBGN option. However, the user 
must refer to the program listing or compilation to determine the meaning of this 
printout . 
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Grid Selection 


The grid selection parameters appear on the third input card and are given by 
DDS, KL, JL, KDS. The number of streamlines including the wall boundaries is given, 
by KL. Good numerical results are obtained with KL = 26. The maximum number of 
streamlines which can be used is KL=100. The number of streamwise stations is 
divided into a coarse grid JL and a fine grid which is JL*KDS. The maximum number of 
streamwise stations in the coarse grid is given by JL=5l. There is no limit to the 
size of the fine grid. Since the numerical method used to integrate the equations of 
motion is conditionally stable, KDS is an important stability parameter. Good results 
are obtained with KDS=5. If the solution is numerically unstable KDS should be 
increased. It should be noted that computing time is directly proportional to KDS.” 
The parameter DDS, distorts the normal coordinate by placing more streamlines near 
the wall. Increasing DDS, moves the streamline closer to the wall. Good results have 
been obtained by setting DDS=10. 

Separation 

Separation is determined when the streamwise component of the wall stress goes 
to zero. Since the equations become numerically unstable as the wall stress 
approaches zero, the wall stress is extrapolated ahead two stations. When the 
extrapolated value of wall stress is less than zero, the calculation is terminated 
and an appropriate message is printed. 
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k.2 Description of Input 


This subsection describes the loading of input data cards for running the 

computer program. Care should be taken in loading the program because the input 

changes depending on the options chosen on the second data card. Multiple cases 

can be run simply by stacking the cases in order. The last case is followed by 

two blank cards.. 


Card No. 

1 : Title 

Card 



Name 

. Col. 

Format 

Comments : 


TITLE 

1-72 ' 

12 Ab* 

Any alpha numeric 

characters .’ 

Card No. 

2 : Option Card 



Name 

Col. 

Format 

Comments : 


I0FT1 

1-2 . 

, 1 

12 

(FLiZfalN Option) 
I0PT=1 The inlet 
data on card 5 « 
I0PT=2 The inlet 
following card 5 . 

flow is computed by specifying the 
flow is read from 2 x KLL data cards 

I0FT2 

2-4 

12 

(F0RCE Option) 



I0FT2=O No blades or struts exist in the duct and 
these cards are not loaded. 

, I0PT2=1 The strut forces are calculated from 
experimental data on cards 2 x KLL cards following 
card 9 . If I0PT1=2 data cards must be separately 
loaded . 

I0PT2=2 The blade force is calculated fromCLandCD 
in Subroutine CASC. 

I0PT2=3 The blade force is calculated from a 3 and Z a 
in Subroutine CASC . 

I0PT2=4 The blade force is calculated from KLL input 
data cards following card 9 with a 3 and Zg given. 

IvPTS 5-6 12 (GDUCT Option) 

I0PT3=1 Calculate a straight annular duct 
I0PT3=2 Head duct coordinates 

I0PT3=3 Calculate a straight wall annular diffuser 
I0PT3=L Calculate NACA curved wall annular diffuser 
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Card No. 2: ■ Option Card (Continued) 


Name 

Col. 

Format 

Comments: = _ ,/^ 

I0PT4 

7-8 

12 

Print solution every I0PT4 stations, i ■ 

I0PT5 

9-10 

12 

Blade force data input (see I0PT2=l) > ' . 

I0PT6 

11-12 

12 

Matching point . - ■ 4 

If not specified, the matching point for Spalding's 
solutions is the second mesh point from the -wall; 
otherwise, it is the I0 Pt 6" I: * 1 point from the wall. 

I0PT7 

13-14 

12 

If I0PT7=O, the gas is assumed to be air. 
If I0PT7^O, input gas parameters. ( • 

Card No. 3 

: Mesh 

Parameters 


Name 

Col. 

Format 

Comments: 

DDS 

1-10 

F10.3 

Mesh distortion parameter (see Section 7*6) 

Good results are obtained 5 < DDS <10. 

■ ; ■ ..v.-.-.-.r 

KL 

11-13 

13 

Number of streamlines including wall 2 5 KL 5 100 

JL 

14-16 

13 

Number of streamwise stations, JLS51 

KDS 

17-19 

13 

Number of steps per streamwise station (see Section 7*6) 

KLL 

20-22 

13 

Number of streamlines of data input (see I0PT1, I0PT2) 

Card No. 4 

: GDUCT 



Name 

Col. 

Format 

\K 

Comments : 


These input cards are read in subroutine GDUCT as 
programmed by the user. 



The following duct geometries (designated as I0PT3=1»2,3j and 4) have been programmed 
Card No. 4; (lgiPT3=l') 


Name 

Col. 

• Format 

Comments : 

RH1 

1-10 

F10.0 

Hub radius (ft) 

RT1 

11-20 

F10.0 

Tip nadius (ft) 

Z1 

21-30 . 

F10.0 

Length (ft) 

TWH 

31-40 

F10.0 

Hub wall temperature (deg R) 

TWT 

41-50 

F10.0 

Tip wall temperature (deg R) 

AMWH 

51-60 

F10.0 

Hub wall bleed (lb/ft^/sec) 

ANWT 

61-70 

F10.0 

Tip wall bleed (lb/ft^/sec) 

Card No. 

4: (l$FT3=2 ) 


Name 

Col. 

Format 

Comments : 

Z1 

1-10 

F10.0 

Duct length (ft) 

Cards Following Card No. 4: 

(I«SPT3=2) 

Name 

Col. 

Format 

Comments : 

R(1,1,J) 

1-80 

8F10.0 

Tip radius (ft) 

R(2,1,J) 

l-8o 

8F10.0 

Hub radius (ft) 

Card No. 

4: (l<zSpT3=R} 


Name 

Col. 

Format 

Comments 

RH1 

1-10 

F10.0 

Hub radius (ft) 

RT1 

10-11 

F10.0 

Tip radius (ft) 

Z1 

11-20 

F10.0 

Duct length (ft) 

ZTHRO 

20-30 

F10.0 

Length to throat (ft) 
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Name 

Col. 

Format 

Comments : 

ANGH 

30-40 

F10.0 

Hub wall .angle (deg) 

ANGT 

30-50 

. F10.0 

Tip wall angle (deg) 

Card No 

. 4: (l<zSPT3=4) 


Name 

Col. 

Format 

Comments : . . _ . 

TOTH 

1-10 

F10.0 

Hub wall temperature, (deg R) 

TWT 

11-20 

F10.0 . 

Tip wall temperature (deg R) 

AMWH 

21-30 

F10.0 

Hub wall bleed (lb/ft 2 /sec) 

AMWT 

31-40 

F10.0 

Tip wall bleed (lb/ft 2 /sec) 

Card No 

. 5:- Inlet 

Flow Distribution 

Name 

Col. 

Format 

Comments : 

AMS1 

1-10 

F10.0 

Nominal inlet Mach number""' - 

ALHL 

11-20 

FI 0.0 

Nominal swirl angle at hub (deg to axis ) 

DSH 

21-30 

F10.0 

Boundary layer displacement thickness on hub wall 

DST 

31-40 

F10.0 

Boundary layer displacement thickness on tip wall 

ANH 

41-50 

F10.0 

( ! ' ; v 

Power law exponent for hub boundary layer 

ANT 

51-60 

F10 .0 

Power law exponent for tip boundary layer 

2 x KLL 

Inlet Flow 

> . h . . . ; - • . ■ 

Cards Following Card 5 fOnlv if Vtiwu?) 

Name 

COl. 

Format 

Comments : 

Y 

1-10 

F10.0 

Fractional distance across duct 

P o 

11-20 

F10.0 

Total pressure (lb/ft 2 abs) 

P 

21-30 

F10.0 

/ ? ’• - -i 

Static pressure (lb/ft abs) 




(ft) 

(ft) 


, i .. 
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a 31-40 F10.0 Swirl angle to axis (deg) 

41-50 F10.0 Total temperature (deg' R) : 

The first KLL cards describe the inlet flow. The 
second KLL cards describe the exit flow. If the exit 
flow is not known , KLL blank data. cards must be used. 
If 8* >0 on Card 5*0> boundary layers are added 
according to card 5 .0. 


Card No. 6: Blade Data (Only if I0PT2/O) 


Name 

Col. 

Format 

Comments: 

ZCL 

1-10 

F10.0 

Axial location of strut centerline (ft) 
(see Fig. (7. 2. 8.1)) 

NB 

11-13 

13 

Number of struts 

ISHARE 

Ik -16 

13 

Blade shape index 


ISHAFE=1 NACA 5 digit series blade 
ISHAFE=2 Swirl vane 

1 


Card No. 7: Blade Row Geometry Hub (Only if I0FT2/O) 


Name 

Col. 

! i 

Format 

Comments : 

RCLH 

1-10 

F10.0 

Hub radius at blade centerline (ft) 
1 

ALPSH 

11-20 

F10.0 

Chord angle to blade face (deg) 

CH0RDH 

21-30 

F10.0 

Chord (ft) 

THICKH 

31-40 

F10.0 • 

Thickness /chord 

! 

PHICH 

41-50 

F10.0 

Camber angle (deg) (see Fig. (7.2. 2.1)) 


Card No. 8: Blade Row Geometry Midpoint (Only if I0PT2/O) 


Name 

Col. 

Format 

c 

Comments : 


RCLM 

1-10 

F10.0 

Midpoint radius 

at blade centerline (ft) 

ALPSM 

11-20 

F10.0 

Chord angle to 

blade face (deg) 

CH0RDM 

21-30 

F10.0 

Chord (ft) 
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Card No. 8 (Continued) 


Name 

Col. 

Format 

Comments : 

THICKM 

31-40 

F10.0 

Thickness /chord 

PHTCM 

41-50 

F10.0 

Camber angle (deg) (see Fig. .(7*2. 2.1)) 

Card No. 

9: Blade Row Geometry Tin (Only if I$FT2/0) 

Name 

Col. 

Format 

Comments : 

RCLT 

1-10 

F10.0 

Tip radius at blade centerline (ft) 

ALPST 

11-20 

F10.0 

Chord angle to blade face (deg) 

CH0RDT 

21-30 

F10.0 

Chord (ft) 

THICKT 

31-40 

F10.0 

Thickness/chord 

PEUCT 

41-50 

F10.0 

Camber angle (deg) (see Fig. (7. 2. 2.1)) 

2 x KLL 

Strut Data Cards Following Card 9 (lfi5FT5&>. I$PT2=1) 

Name 

Col. 

Format 

Comments : 

Y 

1-10 

F10.0 

Fractional distance across duct 

P o 

11-20 

F10.0 

Stagnation pressure (lb/ft 2 abs) 

P 

21-30 

F10.0 

Static pressure (lb/ft abs) 

a 

31-40 

F10.0 

Swirl angle to axis (deg) 

T o 

41-50 

F10.0 

Stagnation temperature (deg R) 


The first KLL cards describe the inlet flow to strut 
row . The second KLL cards describe the exit flow . 
(see Card 2). 
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KLL Blade Data Cards Following Card 9 (I0PT5^O> I0PT2=4) 


Name ' 

'Col. 

Format 

Comments : 

y 

1-10 

F10.0 

Fractional distance across duct 

a 

11-20 

F10.0 

Swirl angle at exit (deg) 

z B 

21-30 

F10.0 

Loss coefficient 

Card No. 

10: Performance Point 

Name 

Col. 

Format 

Comments : 

PKESO 

1-10 

F10.0 

Inlet stagnation pressure (lb/ft 2 abs) 

TEMPO 

11-20 

F10.0 

Inlet stagnation temperature (deg R) 

ACI 

21-30 

F10.0 

Clauser constant (0.016) 

AKI 

31-40 ° 

F10.0 

Von Karman constant' (0.4l) 

AH 

41-50 

F10.0 

Van Driest constant (26.0) 

PRTI 

51-60 

F10.0 

Turbulent Prandtl number 

PR LI 

61-70 

F10.0 

Laminar Prandtl number 

Card No. 

11: Specif 

V Gas (Only if I0PT7=l) 

Name 

Col. 

Format- 

Comments : 

CPR 

1-10 

F10.0 

2 2 

Specific heat at constant pressure (ft /sec ) 

CVR 

11-20 

F10.0 

Specific heat at constant volume (ft 2 /sec 2 ) 

VI SCR 

21-30 

FlO.O 

Viscosity (lb sec/ft 2 ) 
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4.3 Description of Output 

Title Page 

i * 

The output presented on this page is self-explanatory except for the following 
variables •’ 


m i = J g B ? u s dr 

r H 

(4.3.1) 

r Tj 

a, = J g B dr 

r H 

(4.3.2) 

_ | <* r T 

P I = 1717 J 9b^ u s p i dr = PRESI 

H 

(4;3.3) 

r ! 

di = rjr f,\f u s (l ' /2 P u i^ dr = DYNPI 

1 r H 

'• (4.3.4) 

— ! /• r T 1 

= "m; J r 9 B pU s M,dr = MACHI 

r H 

(4.3.5) 

WFL<£ = 32.2 N b M, 

(4.3.6) 

US R = M, /p f /a , 

(4.3.7) 

REY = r r U r p r / yi r 

(4.3.8) 


Wall Conditions Page 

O O 

This page presents a table of Z, Rh » rr> N ,® H , R T , m T , © T ' w ^ ic ^ was calculated 
in Subroutine GDUCT. • 1 
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Strut Data Pages (2) (if I0PT2=l) 


These pages present the strut data that were read into the program with 
I0PT2=1. It is written with the same format as the input (see Section 4.2). 

Blade Force Variables Pages (2) (if I0PT20O) 

a a A A A A n A -A 

These pages present the variables Y,U s ,U^,a, n, ©,M,n 0 ,©o> P across the 
flow passage at the upstream (station 2) and downstream (station 3) control surface. 
These flow variables are used to calculate the blade force. 

Blade Force Page (if I0PT20O) 

This page presents a summary of the blade force calculation Y, G B , H S ,H,£,C L , 
C D ,Z B ,X across the duct. 

Gap Average Flow Page 

This page presents the solution for the flow variables across the duct at 
selected streamwise stations depending on I0PT4. A table of values for Y, U s , U 
a ,11,®, M,n o ,© 0 ,Cp are given where 


C P = III - ri(0,0))/q, 


(4.3.9) 


In addition, the wall values for Z, m,C^,Cj s ,Q are printed, where and 


Cj s are defined by 


= T n0 /q i 


( 4 . 3 . 10 ) 


S =T ns /q i 


( 4 . 3 . 11 ) 


The one -dimensional characteristics of the flow are also given; area ratio a/o 15 
Mach number (isentropic f,low),M| , incompressible and compressible flow pressure 
coefficient CPINC and CPC0MP. 
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Finally, the following quasi -two-dimensional boundary layer information is 
given : 
t 


THET = Q*/r f 

H12 = 8 */B* 

REYTH = 'co u co 9 " / '‘® 

cf = V^ao 

REC = R c 

STAN = St 


, Momentum thickness (dimensionless) 

, Shape factor (dimensionless) 

, Reynolds number based on momentum, thickness (dimensionless)- 
, Wall friction coefficient (dimensionless) 

, Recovery factor (dimensionless) 

, Stanton number (dimensionless) 


These variables are based on boundary layer thickness defined by the maximum 
velocity in the duct and, therefore, may not be meaningful if swirl, flow distortion, 
or large normal pressure gradients exist. 

IDBGgi Pages . 

' Intermediate printouts which were used to debug the program may be called by 
setting the debug options IDBG0=1. IDBG0 may be specified on the option input card. 
The user should refer to the program listing in each subroutine to determine the 
printout variables . 


IDBG0 

Subroutine 

Printout 

IDBG1 

SPLDC 

Tables of U + (Y + ) , H * ( Y + ) , H + (Y + ) ,.W + ( Y + ) 

IDBG2 

SPLDC 

Table of C N (m + ) , C N (P + ) 

IDBG3 

FL0WIN 

Flow input variables and parameters 

IDBG4 

STRESS 

Table of G(K) and parameters 

IDBG5 

S0LVE 

Tables of DF(K), DHF(k), DPF(k) 

IDBG6 

C00R 

Tables of Q(k), RH(j), RT(J) 

IDBG7 

F0RCE 

Table of FF(k) 
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• ■ -4.4 Description of Failure Modes . > 

The computer program can diagnose the cause of certain failure modes and a 
printed message of the following form is given. 

**DIAGNOSTIC NO . <t> FOR ANNULAR DIFFUSER DECK** 

The number 0 identifies the type: of failure from the list below. 

(1) AMPLUS OR PPLUS OUTSIDE RANGE OF TABLE 

This failure occurs in Subroutine CINTP. In attempting to interpolate a value 
of C(m*) or C(p + )from the table, a value of m* or p + is encountered which is outside 
the tabulated values. Calculation is then terminated. This failure may occur for 
small values of U* as the flow approaches separation. 

r • 

(2) ITERATION IN FUNCTION USTAR FAILED TO CONVERGE 

This failure occurs in Function USTAR and means that Newton's iteration does not 
converge. Since U + CY + ) is a well behaved function with no saddle points or 
singularities except a logrithmic singularity at Y* = 0 , Newton's method should 
converge if the input is correct. The possibility of extrapolating to a negative 
value of U* has been anticipated and corrective measures applied. Failure, however, 
may occur when the numerical integration is numerically unstable, or if Y* < 10 . 

(3) AMPLUS CANNOT BE LESS THAN -0.06 

This failure mode occurs in Subroutine SPLD because asymptotic solutions do not 
exist for m 0.06 . Calculation is terminated. 

(4) ADD STRAIGHT ANNULAR CHANNEL INLET 

I 

This failure mode occurs in Subroutine C00R and represents a warning that the 
inlet flow must not have a normal pressure gradient. For better results add a 
straight annular channel inlet. The calculation, however, is not terminated by this 
failure mode . 

(5) WALL CURVATURE TOO LARGE AT STATION (i) 

This failure mode occurs in Subroutine C00R when the wall curvature is too large 
and the calculation is terminated. This failure frequently occurs when the duct 
coordinates are discontinuous. Check the calculation in Subroutine GDUCT. 
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(6) SOLUTION UNSTABLE AT S = (XXX.) 

This failure mode occurs in Subroutine S0LVE when the streamwise velocity is 
negative resulting from numerical instability. The calculation is then terminated. 
Decreasing the step size (increasing KDS) may improve stability. 

(7) BLADE DATA IN ERROR IN FORCE INTERPOLATION ROUTINE 

This failure mode occurs in Subroutine ,F0RCE from failure of the interpolation 
procedure. This failure can only occur if the input data is wrong. 

(8) ITERATION OF BACK PRESSURE CALCULATION FAILS TO .CONVERGE 

This failure mode occurs in Subroutine FINVIS when the iteration scheme required 
to compute the back pressure from radial equilibrium of the downstream flow fails to 
converge. This failure frequently occurs when the loss coefficient is excessively 
high. Calculation is terminated. 

(9) YPLUS LESS THAN 10. IN FUNCTION FUPLUS 

This failure mode occurs in Function FUPLUS because the inner layer analytical 
solutions are not valid for Y + < 10. This failure sets U +B 0 and, therefore, the 
calculation must be terminated. The problem may be corrected by increasing Y* at 
the matching point, by increasing the Reynolds number, or increasing the mesh spacing 
near the wall. 

(10) YPLUS AT MATCHING POINT OUTSIDE BOUNDS ,, 

For best results, Y + at the matching point must lie between the limits 40 < 

Y*< 500. When the calculation is outside these bounds, Y is printed out but the 
calculation is not terminated. 

(11) PROGRAM ASSUMES INLET FLOW HAS CURVATURE 

The analysis assumed that the inlet flow duct has no curvature, hence normal 
pressure gradients are produced only by swirl. Calculation, however, is not terminated. 

(12) CRITICAL MACH NUMBER INSTABILITY 

This failure mode occurs when the critical Mach number is reached too close to 
a mesh point. It can be corrected by changing KDS to move the mesh point. 


(13) ISHAFE=(l ) AND I0PT2=(j) ARE NOT CONSISTENT 

The input options are not consistent. See Section 4.2. 
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6.0 LIST OF SYMBOLS 


a Area (ft 2 ) 

2 

A Area, a/r r (dimensionless) 

A* Van Driest constant (26.0) 

A* Critical area ratio (dimensionless) 

b Chord (ft) 

B Chord, b/r r (dimensionless) 

C Speed of sound (ft/sec) 

Cp Specific heat pressure (ft 2 /sec /deg R) 

C v Specific heat volume (ft 2 /sec 2 /deg R) 

2 

Cp ' Drag coefficient, 2D/(p2 u 2 (dimensionless) 

2 

C L Lift coefficient, ZL/ipz^z b) (dimensionless) 

.0 Drag/span (lb/ft) 

e ns ,Streamwise strain (l/sec) 

e n< £ Tangential strain (l/sec) 

f Force/area (lb/ft 2 ) 

f Force /span (lb/ft) 

g Gap between walls (ft) 

G Gap between walls, g/r r (dimensionless) 

9e Gap between chord lines (ft) 

Gg Gap between chord lines (dimensionless) 
h Enthalpy (ft 2 /sec 2 ) 


6l 



Universal stagnation enthalpy, (how "ho)pw U*/q w.(dimensionless ).. . 
Universal adiabatic stagnation enthalpy, (h 0 A W -h 0A ) /(U*) • (dimensionless) 
Entropy (ft 2 /sec 2 /deg R) .. ... 

Lift/span (lb/ft) . .. . - 

Mass flow (slugs/sec) 

Mass flow, mANgiy 2 ^ U r ) (dimensionless) __ ... .... 

Mass flow/area (slugs/ft 2 /sec ) . 

O 

Universal mass flow parameter, m w /{p w U*). (dimensionless) 

Mass flow/area, m/(p r u r ) (dimensionless) 

Mach number, U/c (dimensionless) ... ..... 

Streamwise Mach number (dimensionless) • 

Normal coordinate (dimensionless) . - „ 

Normal coordinate, n/(r r V r ) (dimensionless.) 

Number of struts (dimensionless) 



Heat flux, q /{p r U r C p T f ) (dimensionless) 
Universal heat flux, q/q w (dimensionless) 
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Universal heat, flux (adiabatic), q}(p w ll* 3 ) (dimensionless) 
r Radius (ft) 

R c Recovery factor, Eq. (3*2.35) (dimensionless) 

R Radius, r/r r (dimensionless) 

(R Gas constant (ft 2 /sec 2 /deg R) 

S Streamwise coordinate (dimensionless) 

S Streamwise coordinate, s/(r r V r ) (dimensionless) 

St Stanton number (dimensionless) 

t Blade thickness (ft) - 

T Temperature (deg R) 

T + Universal temperature, CpT/u* z (dimensionless) 

u s Streamwise velocity (ft/sec) 

u n Normal velocity (ft/sec) 

Tangential velocity (ft/sec) 
u Magnitude of velocity (ft/sec) 

Ug Blade velocity (ft/sec) 

u* Friction velocity, -/r w /p w (ft/sec) 

U s Streamwise velocity, U s Aj r (dimensionless) 

Un Normal velocity, U n Aj r (dimensionless) 

U<£ Tangential velocity, U(£/Ur (dimensionless) 

U Magnitude of velocity, U/U r (dimensionless) 

U B Blade velocity, U B /U r (dimensionless) 
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Universal velocity, U/U* (dimensionless) 

Friction velocity, U*/U r (dimensionless) 

Metric scale coefficient (dimensionless) '■ 

Metric scale coefficient, v/v r (dimensionless) 

Stream function inner layer (dimensionless) 

Volume (ft3) 

Distance along streamline (ft) 

Distance normal to wall (ft) 

Axial distance (ft) 

Distance along streamline, x/r r (dimensionless) 

Distance normal to wall, y/r r (dimensionless) 

Universal distance from wall, y p w U*/f± w (dimensionless) 
Axial distance, z / fr (dimensionless) 

Loss coefficient (dimensionless) 

Swirl angle to axis (deg) 

Chord angle to axis (deg) 

Ratio of specific heats, Cp/C v (dimensionless) 

Boundary layer thickness (ft) 

Displacement thickness (ft) 

Boundary layer thickness, 8/r r (dimensionless) 

Displacement thickness, SVr r (dimensionless) 

Streamwise strain, r r e ns /u r (dimensionless) 

Tangential strain, r r e n<jb/u r (dimensionless) 



V 


Transformed normal coordinate (dimensionless) 

H Blade/force area r r f/[p T U r ) (dimensionless) 

8 Angle of streamline to axis (deg) 

9* Momentum thickness (ft) 

© Temperature ratio, T/T r (dimensionless) 

©* Momentum thickness, 9*/J f (dimensionless) 

I Entropy, (X-I r )/(R (dimensionless) 

k Von Karman constant (0.4l) 

X Thermal conductivity (lb/sec/deg R) 

p. Viscosity (slugs/ft/sec) 

0 Blade force/span, $ /(r r p r u r ) (dimensionless) 

II Pressure ratio, p/p r (dimensionless) 

p Density (slugs/ft^) 

p s i p^ Radius of curvature (ft) 

P Density ratio, p / p f (dimensionless) 

P s ,P n Radius' of curvature (dimensionless) 

cr Solidity, b/g B (dimensionless) 

ns Streamwise stress , T n s/(/3r u r ) (dimensionless) 

Tangential stress, r rs( p/{p r [i x z ) (dimensionless) 

r ns Streamwise stress (lb/ft 2 ) 

T n ^ Tangential stress (lb/ft^) 

T * Stress, r /t w (dimensionless) 
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<jf> Tangential coordinate (radians) 

<p c Camber angle (deg) 

Blade dissipation function (lb/sec/ft 2 ) 

<f> B Blade dissipation function (dimensionless) 

X Clauser constant (0.016) (dimensionless) 

X Normal coordinate transform, d^/dn (dimensionless) 

Stream function (slugs/ft) 

'f' Stream function (dimensionless). 

Subscripts 

0 Stagnation conditions 

1 Inlet conditions 

2 - Upstream of strut 

3 Downstream of strut 

A Adiabatic 

E Effective turbulent 

H Hub conditions 

1 Incompressible conditions 

M Midspan conditions 

r Reference conditions* 


*Reference conditions are based on standard sea level atmosphere conditions for all 
thermodynamic quantities. The reference radius r r is the inlet outer radius, and 
the velocity is the mean inlet velocity. 
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T Tip conditions 

w Wall conditions 

CD Free stream or edge of boundary layer 

Superscripts 

— Mean or average quantity 
A Variables for blade force calculation 

/ Deviation from mean quantity 
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7.0 APPENDIX A: DETAILED PROGRAM DESCRIPTION 


A detailed description of the computer program is given in this section. The 
subroutines and external functions are described indi'vidually. A list of F0RTRAN 
variables and DIMENSION FORTRAN variables in blank common are given as well as a 
source listing and a set of functional flow charts. 

The subroutines and external functions are all described with the same format 
using the name of the subroutine with Its argument list given as a title. A list of 
options and F0RTRAN symbols' used only in the named subroutine are; then given. Any 
special or additional theory used in the subroutine is presented but well known 
numerical methods are not described. Figures and Flow Charts appear with the 
subroutine. 

Subroutine GDUCT, CASC, and FTHIK, which generate the duct shape, and calculate 
the strut lift and drag coefficient, and blade thickness distribution, respectively, 
have special status since they must be written by the user . Sample calculations are 
given, and the user should examine carefully the flow charts for these subroutines. 
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7.1 Main Rrogram 


Object Control input/output . 

Options All options described In Section 4.2. 

List of Symbols All variables In blank C0MM0N. 

Theory 

The MAIN" program serves as an executive program to read input according to 
selected options and write the appropriate output. In addition, it computes or sets 
fixed data which is used for any case that may be calculated such as calling 
Subroutine SPLD which computes the constants for the inner layer solution. Finally, 
the MAIN program calls Subroutine S0LVE which does the actual computation. The flow 
chart is shown in Fig. 7.1.1. 



/ READ 
7 OPTIONS 


/read/write/ 

/ MESH h 
3 ARAMETERs/ 


WRITE 

title 

. OPTIONS 



/ READ DUCT I 
GEOMETRY / 


r PROGRAMMED BY USER 
IN SUBROUTINE GDUCT 


/READ/WRITEy 

/inlet flow/. 

PARAMETERS/ 


/ READ INLET/ 


/READ/WRITE 
'blade data/ 


/ BLADE 
/ DAT A = 
HNLET F LOW y 
DATA / 


BLADE FORCE 
COMPUTED FROM 

INLET/EXIT 

LY,P n .P^.T n 


/READ BLADE/ 
/FORCE DATA/ 


' read/write 

PERFORMANCE, 
POINT / 


/WRITE DUCT y 
7 SHAPE / 
GEOMETRY / 


> WRITE INLET, 
BLADE FORCE 
Y,P 0 ,Pa.T 0 / 


/ SET UP/ 
'CALCULATE, 
FLOW / 


FIG. 7.1.1rFL0W CHART MAIN PROGRAM. 


70 





7.2 List and Description of Subroutines 



Name 

Ob.lect 

7.2.1 

BL2D 

Compute Boundary Layer Parameters 

7.2.2 

- CASC 

Compute Strut Performance 

7.2.3 

CINTP 

Interpolate Table of Asymptotic Constants 

7.2.1+ 

C00R 

Compute Coordinate Functions 

7.2.5 

FINVTS 

Calculate Strut Back Pressure 

7.2.6 

FL(ZfalN 

Setup Inlet Flow • 

7.2.7 

F0RCE 

Calculate Strut Force 

7.2.8 

GBLADE 

Calculate Blade Geometry 

7.2.9 

GDUCT 

Compute Duct Shape 

7.2.10 

SLETE' 

Find Strut Control Surfaces 

7.2.11 

S0LVE 

Integrate Equations of Motion 

7.2.12 

SPLDC 

Compute Constants for Inner Layer Solution 

7.2.13 

STRESS 

Compute Stress and Heat Flux 

7.2.14 

VARB 

Compute Blade Force Variables 
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7.2.1 Subroutine BL2D 


Ob.j ect Compute Boundary Layer Parameters 
Options 


11=1 Calculate Hub Boundary Layer 


111= 

-2 Calculate Tip Boundary Layer 
=1 Curve Fit Y = PU /(PcoLIqq) 

=2 Curve Fit Y =U/Uco 

List of 

Symbols 

• 

AMF 

= M co 

5 . • ‘ . 

, Free stream Mach number (dimensionless) 

AN1,AN2 

= n | ,n 2 

, Velocity Profile Power Law Exponent (dimensionle 

BL0CK 

= B c = [a£ R h + AVR t ]/A 

, Blockage factor (dimensionless) 

CF1 

= Pqo u cd^ 

, Wall Friction Coefficient (dimensionless) 

DSH1 

' a h=^ Ah 1 | -|d^]'<>' 

=av/ o a V|^w 

0 r QOQO 

= A7A , 

, Displacement thickness hub (dimensionless) 

DST1 

D1 

, Displacement thickness tip (dimensionless) 

, Displacement thickness ratio (dimensionless) 

D2 

= ©7a. 

, Momentum thickness ratio (dimensionless) 

H12 

= H l2 = A*/©* 

, Shape factor (dimensionless) 

QF 

■ “cO = 1/2 P 00 U ® 

, Free -stream dynamic pressure (dimensionless) 

QW 

= Q oj 

, Wall heat flux (dimensionless) 

RCALC 

= R c -[®w©CO' 1 ] /(^ 2 McO; 

), Recovery factor (dimensionless) 

BEYTH 

= R e0= PooHe ©V;® 

, Reynolds number (dimensionless) 

RF 

-Pco 

, Free jstream, density (dimensionless) 

STAN 

= Q w f [ffo u a/® w - ®W] 

, Stanton number (dimensionless) 

TF 

ii 

© 

8 

, Free -stream temperature (dimensionless) 
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THET 

= ©* = / 

,A Pu s r 
5 PUoo 1 

u s I 

U CD J 

l dY 

, Momentum thickness (dimensionless) 

TW 

n 

© 

< 




, Wall temperature (dimensionless) 

TWA 

ii 

r, + P |/3 J^j 

L RL 2 

Lm 2 

M oo 

© 

8 

, Adiabatic wall temperature (dimensionless) 

UF 

= Uco 




Free -stream velocity (dimensionless) ; ' 

XB 

= X 




, Average X (dimensionless) 

YB 

= Y 




, Average Y (dimensionless) 

Theory 







When the boundary layers are very thin, it is difficult to integrate the flow 
■variables accurately. Therefore, a power law is fitted to the solution of the form 

l/ni 

PUAPcoUcq) = (Y/A,) (7.2. 1.1) 


and 

■i 

U/Uqo = (Y/A,) l/n2 -(7.2. 1.2) 


which with logrithmic differentiation become 

ln(y/A,) = n, In [lPU)/(Poo Uoo)] (7. 2. 1.3) 


or 

' ! ln(y /A 2 )= n 2 ln [u/U^] 

• i r •>. 


1 (7. 2.1.4) 
. / 
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Letting 


X = ln(Y/A,) or ln(Y/A 2 ) (The choice dependent upon III) 


and 


'n( P V P ooU®]" ,l "[ U s /U co] 

a least square fit 


Y=Y + b(X — X) 


is used to determine n ( ,n 2 ,A| , A 2 • Finally , 


A*/A ( = 1/(1 + n ( ) 


and 



( 7 . 2 .' 1 . 5 ) 


. * ,■ ' * <■ 
iV » - 

’( 7 . 2 . 1 . 6 ) ‘ 


( 7 . 2 . 1 . 7 ) 


( 7 . 2 . 1 . 8 ) 


( 7 . 2 . 1 . 9 ) 



: 7.2.2 Subroutine CASC (Argument List) 

Ob.iect Compute Strut Performance 

' * ■ * * . . - . 

Options 

ISHAPE=1 NACA 5 digit Series Airfoil 
=2 Flow Turning Vane 

List of Symbols (Argument List) 


ALPH3 

=-a 5 

j Exit Flow angle (deg) 

ALPS 

= a s 

, Chord angle to axis (deg) 

ALP1 

A 

* a 2 

, Inlet flow angle (deg) 

AMACH 

= Cig 

, Inlet Mach number (dimensionless) 

BCH0RD 

» B 

, Chord (dimensionles^ 

CD 

II 

o 

O 

, Drag coefficient' (dimensionless ) 

CL 

■ = c L 

, Lift coefficient (dimensionless ) 

PHIC 

3 4>c 

, Camber angle (deg ) 

S0LD 

= a- « B/G b 

, Solidity (dimensionless* 

TM 

= 

, Blade thickness (dimensionles^ 

ZLjzSSS • 

= 2 b 

, Loss; coefficient ' (dimensionless ) ' 

* ’ 1 * ’ ■' ‘ ‘ 

Theory 




This subroutine computes the blade or strut section performance. Since many 
strut shapes are possible (Refs. 10 and 11), only two shall be described here for 
illustrative purposes. 

ISHAHXL 

The subsonic performance of isolated airfoils is given in Ref. 11. Selecting 
the NACA 5 digit uncambered airfoils and applying the Prandtl-Glauert rule 
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dc. 


da /o 


t 


M 


= 0.109 - 0.0256 


B V \- rOl | 


( 7 . 2 . 2 . 1 ) 


e, = 


dc L \ (a 2 -a s ) 


L V 


( 7 . 2 . 2 . 2 ) 


C 0 = 0.0075 - 0.00382 C|f 


(7- 2. 2. 3) 


ISHAEE=2 

For turning vanes of high solidity, it is assumed that the exit flow angle is 
equal to the metal angle and that the loss coefficient is given by the stator cascade 
losses of Ref. 12. Then 


A 

a 


3 


a s~4>c /Z 


(7. 2. 2. 4) 


D 



A 

cos a 2 
cos a 3 


] 


+ 


A 

cos a 2 

2 cr 



(7. 2.2.5) 


4 28 

$ D = 0.005 + 0.005 (D/4) ' 


(7- 2. 2. 6) 


Z B = <D 


2 cr f COS 02 ] 2 


D COS 3 3 l COS 


]' 


(7. 2. 2. 7) 
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<2 


' The' input to this subroutine comes in the argument list and is given by a 2 , a s , 
2 , B, <f> c , t M , c ISHAPE, I0PT2, N0PT5 • The output is also in the argument list and 
is either C L , C D , or a 3 , Z B . It is noted that I0PT2 and ISHAE2 must be compatable. 
The sign convention for determining this force is shown in Fig. 7.2. 2.1 and the flow 
chart for this subroutine is given in Fig. 7*2. 2. 2. 



V 


z 

1. e DECELERATING FLOW 


Z 

1. d ACCELERATING FLOW 



FIG. 7.2-2.1 -SIGN CONVENTION BLADE FORCES. 
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FIG. 7.2.2.2-FLOW CHART FOR SUBROUTINE CASC. 
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7-2.3 Subroutine CINTP 

. 4 

Object Interpolate Table of Asymptotic Constants 

y • 

Options 

IF(l < LO < LM-2 ) Calculate Constants 
IF( LO^l. OR. L0=»LM-2) Stop Calculations 


List of Symbols 
AMPLUO = mo* 

LO 

PLUSO = p 0 + 

X 

Theory 


, First value of interpolation for m + (dimensionless) 

, Index of first value of interpolation 
, First value of interpolation for p + (dimensionless) 
, Fractional distance (dimensionless) 

i 


Newton's forward differencing interpolation formula is used since it involves 
only tabulated values. Let 


x =(m*-mJ)/Am t 


(7-2. 3-1) 


= (P + -P*)/A P + 


(7. 2. 3. 2) 


Then for three successive tabulated points ? for any f^ulated function f 


} = fo + (frJo) * + ^ 2 “ 2/ ( + / 0 )x(x-l)/2 


(7. 2. 3. 3) 


with the restiction for best accuracy that 


X < I 


(7.2.3 .4) 
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7.2.4 Subroutine C00R(KSS, JSS) 
Ob.iect Compute Coordinate Functions 
Options 

N0PT3=1 '-Compute and Store Coordinates 
=2 Compute Coordinate Functions 


List 

of Symbols 


A1 

= R t (J+I)-R t (J“I) 


A2 

= R t (J+I)-2R t {J) + R t (J-I) 


A5 

= Y h +Y T 

, Sum of metric scale coefficients (dimensionless) 

a6 

= V V T 

, Ratio of metric scale coefficients (dimensionless) 

B1 

= R h (J + I)-R h (J-I) 


B2 

= R H (J+l)-2R H (J)+R H (J + n 


C 

= c 

, Constant in x transformation (dimensionless) 

Cl 


• « 

, Mean wall slope at J-l (dimensionless) l 

C2 

= Z T*~ Z H 

, at J-l (dimensionless) 

03 „ 

= r t“ r h 

, at J-l (dimensionless) 

c4 


, Slope difference at J-l (dimensionless) 

D1 


, at J+l (dimensionless) 

D2 

X 

M 

1 

»- 

M 

II 

, at J+l (dimensionless) 

D3 

II 

-T 

, at J+l (dimensionless) 

d4 

= / dStt. + / 2 + hzh 
\ dZ dZ ) /c R t -R h 

, Slope difference at J+l (dimensionless) 

DZ 

= AZ 

, Step size (dimensionless) 

RH1 

= r h 

, Wall radius at s hub (dimensionless) 
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RH2 

= dR H /dz. 

, Wall slope at s hub (dimensionless ) .. , 

RH3 

= d 2 R H /dz 2 

, Wall second derivative at s hub (dimensionless) 

RT1 

= R r 

, Wall radius at s tip (dimensionless) 

RTS 

= dR T /dz 

, Wall slope at s tip (dimensionless) 

RT3 

= d 2 R T /dz 2 

, Wall second derivative at s tip (dimensionless ) 

YT 

= y t 

, Duct height at S (dimensionless) 

ZH 

11 

M 

I 

, Axial location of hub wall at s (dimensionless) 

ZT 

= z T 

, Axial location of tip wall at s (dimensionless) 


Theory 


An approximate potential flow solution can be constructed for ducts where the 
meridonal plane curvature and divergence of the duct is not too large (see Section 3.1). 
This solution for the coordinates must be known in terms of n and s rather than R 
and Z since these are the independent variables. 


Through each mesh point on the hub contour a straight line is passed, through 
the tip contour such that its slope satisfies the condition 


^ r m _ | r dR H dR T I _ Z T Z H 

dz 2 i dz + d z J = r t -r h 


(7.2.4. l) 


where R T and Z T are the intersections as shown in Fig. 7*2. 4.1. This straight line 
then approximates the potential line. Since the metric scale coefficient is the 
potential flow velocity, it follows that 



V H + V T 
2 



( 7 . 2 . 4 . 2 ) 


V H _ ^ X T 

V T ' . AX h 


(7. 2. 4. 3) 
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The streamwise coordinate, from Eq. (3.1.8) is 


s = / vdx 

*'o 


(7.2 A. 4) 


The normal coordinate is transformed in the following manner so as to place more 
streamlines near the wall 


n = + tan' 1 [2 tan -y (2 - )] 


(7. 2. 4. 5 ) 


The constant C is a mesh distortion parameter which can be varied between the limits 


0 < C < 77- 


( 7 . 2 . 4 . 6 ) 


The coordinate functions are constructed as follows 


R(n,s) = R H (s) + [r t (s) - R h (s)] n 


(7. 2. 4. 7 ) 


and 


Z(n,s) = Z T (s) + [z T (s)-Z H (s)]n 


( 7 . 2 . 4 . 8 ) 


For the special case of no centerbody, the following transformation is used 


n = tan' 1 [(tan c ) 77] /c 


(7. 2.4. 9 ) 
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c 


(7.2.4.10) 


tirj 

dri 


d Z y _ 2C 2 sin (cn) _ 2 C sin 

dn 2 tone cos ? (cn) cos 


(7-2.4.11) 
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T 


z 


FIG. 7.2.UrCONSTRUCTION OF COORDINATES. 
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7.2.5 Subroutine FINVIS(L0PT) 


Object Calculate Strut Back Pressure 


Options 


L0PT=1 Input Upstream Flow 

=2 Calculate n 3 (0) from °3> Z B 
=3 Calculate Downstream Flow 
=4 Calculate ^ ^ from C L , C D 


List of Symbols 

DFXIC 

= 

An 

EPS 

= 

€ 

KIMH,KIMT 

= 


PI HI 

= 

Am 

ny<0) 

PIH2 

= 

n l ‘ +l \o) 

PIMIN 

= 

n M 

PSIHH 

= 

o 

ro 

PSIHT 

= 

A 

V'a (') 

PSIH1 

= 

A(l) 

PSIH2 

= 


Theory 




, Increment in back pressure (dimensionless) 

, Convergence criteria (dimensionless) 

, Index for end points; hub, tip (dimensionless) 

, Iterative guess for back pressure (dimansionless) 
, Iterative guess for back pressure (dimensionless) 
, Minimum total pressure (dimensionless) 

, Stream function at hub (dimensionless) 

, Stream function at tip (dimensionless) 

, Iterative solution of $ at tip (dimensionless) 

, Iterative solution of $ at tip (dimensionless) 


This subroutine calculates the upstream and downstream flow variables needed to 
determine the blade force and dissipation function, as shown in the flow chart 
(Fig. 7. 2. 5.1) 

Lg>PT=l 

3 

The upstream flow is assumed to be equal to the computed flow at the upstream 
control surface. In the annulus boundary, layer , the variables are set equal to their 
value at the boundary layer displacement height (see Section 3.3). 
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L0PT=2 


A 

The downstream back pressure n 3 (0) , is determined by integrating Eq. (3-3 .30) 

and Eq. (3-3.31) using a third order Runge-Kutta integration. The initial value 

n ( 0 ) must be determined iteratively (see flow chart, Fig. 7. 2. 5-2). A pressure 

increment Afl is estimated from the minimum total pressure n as follows 

M 

An = n M /50 ■ (7. 2 . 5.1) 


Then guesses for the hub wall pressure are obtained from 


A A 

n 3 to) * n M +(An)dTER) 


(7.2. 5. 2) 


Successive guesses are made for 11,(0) with! An increments, until the tip boundary 
a J A(i-D A(i-2) 


condition ((/ (I) is bounded by successive guesses y, (|) ( \^/ (|) • Then using 


Newton's method 


A o r 3 

n,ioj = n,(o) + 


■A A i 

+M) ~ ( I) 


^ ( 2 ) (l)-^ 0 (l) 

■3 o 


[n l 3«» - n ( 3(oi] 


(7. 2. 5. 3) 


Convergence occurs when 


njlo) - fi^o) ) / n ( f(6) 


< e 


(7. 2. 5. 4) 


L0PT=3 ■ 

A ' A 

A When n 3 (T 7 ) is known, then the downstream flow can be determined since Uq^tj 
and © 3 (t)) is known (Eqs. (7-2. 7.1) through (7-2. 7. 7)). 


LgSPT=4 


This option integrates Eq. (3*3.28) and Eq. (3*3.29) using a third order 1 
Runge-Kutta integration. Subroutine VARB is used’ to ; get the integrands for the 
given flow variables . An iterative procedure similar to L0PT=2 option is used to 
determine FI 3 ( 0 ) . 
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FIG. 7.2.5.1.-FLOW CHART FOR SUBROUTING FINVIS. 
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BACK PRESSURE ITERATION 



FIG. 7.2.5.2-FLOW CHART FOR BACK PRESSURE ITERATION LOOP 
IN SUBROUTINE FINUIS. 
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7-2.6 Subroutine FL$tfIN 

Object Set Up Inlet Flow 
Options 

I0FT1=1 Compute Inlet Flow 

IF (T 0 = 0) Ffe = P r and T 0 = T r 

I0PT1=2 Read Inlet Flow 

IF ( 8 >0) Add boundary, layers specified;, 
IF (T >10.) Normalize with Tr and Pr 

List of Symbols 

BINP(l j J,K) , Interpolated from B:ENPUT(lH,J,L) 

Theory ■ I0PT1=1 , 


For this option the free -stream flow is assumed to be isentropic with a constant 
free -stream Mach number M s and wall boundary layers defined by power law velocity 
profiles. Since the swirling flow must be in- radial equilibrium, Eq. (3.2.39) must 
be satisfied together with the isentropic flow relations, Eqs. (3.2.51) and (3.2.52). 
Neglecting curvature in the meridonal plane 


an . y_ 

an " R : 



(7. 2. 6.1) 


where 

'■ 2 2 ' 2 
= M - M s 


(7. 2. 6. 2) 


r- 1 





(7. 2. 6. 3) 


/ 
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Equation ( 7 . 2 . 6 . 1 ) can be integrated with 


and 


m( 0) = M s /cos (a H ) 


( 7 . 2 . 6 . 4 ) 



(7. 2. 6. 5) 


as initial conditions using Runge-Kutta method. 

For a given displacement thickness and velocity profile power law, wall boundary 
layers can be added, assuming collateral boundary layers, such that a is unchanged. 
Then 


A= ( i + n) A* 


( 7 . 2 . 6 . 6 ) 


_LL = /J_Y /n 2 
Uoo \ A / 


(7.2. 6.7) 


and 


©_ 

@00 


I + -n/Prl ^~2~ M C 0 


HI]] 


U Yh + ®w;®AW , , 


© 


00 




— 1 
Uoo J 


( 7 . 2 . 6 . 8 ) 


Finally, the inlet mass flow and reference velocity are determined as follows 



(7.2. 6. 9) 


w = Qp r U r Q r A (7.2.6.10) 
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Theory - ldPTl-2 


For this option, the input flow is calculated from experimental input data. 

The- input variables selected are spanwise location, total pressure, static pressure, 
flow angle, and total temperature, since these are the primary measured variables. 

A simple linear interpolation is used so that for any variable , 

<£(Yj[i7)) = <£(Y|) + [<£(Y 2 )-<£(Y|)][ Y(t7)“ Y|] (7.2.6.11) 


The flow variables are calculated from Eqs. (7. 2. 7.1) through (7*2. 6. 7). 

If ( T 0 >lO), it is assumed that pressure and temperature are given in PSF and 
deg R, respectively, and the flow is normalized accordingly. If 8* is given, it is 
assumed that boundary layers should be added accordingly to the velocity profile 
power law above. Finally, the weight flow and reference velocity are determined 
fromEqs. (7. 2. 6. 9) and (7.2.6.10). 
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FIG. 7.2.6.1-FLOW CHART FOR SUBROUTINE FL0WIN. 
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7.2.7 Subroutine F0RCE(L0PT) 


Object Compute Blade Force 
Options 

L0PT=1 Compute Blade Force/Span 
=2 Compute Local Blade Force 

I0PT2=1 Use Experimental Data 

=2 Use Isolated Airfoil Predicted C L and C D 
=3 Use Cascade Predicted Zg and < 2 % 

=4 Use Input Zgand a 2 

List of Symbols 

AINPT(I,J,K), Interpolated From AINPUT( I +1,J,L) 

Theory LdPT=l 

The theory is given in Section 3*3 and a flew chart in Fig. 7. 2. 7.2. 
IgfeT2=l 


This option computes the strut force from experimental blade element data taken 
just upstream and downstream of the strut. If y, n o ,.n,© 0 , Cl are measured at the 
two locations, we have , 



(7.2. 7.1) 



A 

P = 


A A 

n/@ 


(7. 2. 7. 3) 


9 k 




(7.2 .7.4) 



or 

A 

Inn 


(7. 2. 7. 5) 


A A a 

u s = u cos a 


(7.2. 7. 6) 


= u sin a' " (7. 2. 7. 7) 


The strut ’force is computed from Eqs . (3.3.18) and (3.3.19) and the losses from 
Eqs. (3.3.22) and (3*3.23). Since the experimental measurements do not necessarily 
occur on streamlines used for computing, a linear interpolation is used to obtain 
AINPT from AINPUT. > ' ' ‘ 

Ig$PT2=2 

If the solidity is small, then the blade or strut force can be calculated from 
Isolated airfoil data in which C L and C D are presented for various blade shapes and 
incidence angles. When this is possible we have, from Eqs. (3.3.25) and (3.3.26) 

S CT = ~ G B^z^ 2 a ' [ _ c l sin a 2 + C D cos a 2 ] ( 7 . 2 . 7 . 8 ) 


H^= “GgPgUg o- '[ -’C L cos a 2 + C D sin a 2 ] (7. 2. 7. 9 ) 
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I(Z)PT2=3 

If the solidity is large, then the blade forces must be obtained from cascade 
data because the blade or struts no longer act as Isolated airfoils. Since these 
correlations are usually given for Z B and a 2 > Eqs. (3.3.30) and (3*3*31) are used to 
satisfy mass flow and radial equilibrium' using Subroutine FINVIS. The force is 
calculated fromEqs. (3.3* 18 ) and (3.3.19)* 

I(Z)PT2=4 

This option is the same as I0PT2=3 except that Z B and a 2 are obtained from 
experimental data input. 


LriPT=2 


This option provides the local blade force once the blade force per unit span 
is calculated. Obtaining the blade force per unit span from Eqs. (3*3.20) and 
(3.3*21), the local blade force is determined as follows. 


H 


s 



(7.2.7*10) 


a 4> 

H< * = ~g b (z te -z le ) 


(7*2.7*11) 
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FIG. 7.2.7.1.-FLOW CHART FOR SUBROUTINE F0RCE. 

























7.2.8 Subroutine GBLADE(KK,L0PT) 


Ob.lect Compute Blade Geometry 
Options 

L0PT=1 Compute Blade Centerline Properties 
=2 Compute Local Properties 

List of Symbols 


RR 

y cl 

> Fractional distance along blade centerline (dimensionless) 

THICKN = 


j Maximum blade thickness (dimensionless) 

ZBAR * 

2 

, Fractional distance along- chord line (dimensionless) 

KK 


Index of spanwise station 

Theory 




If three points are given for any blade parameter along the blade centerline, a 
parabolic fit is used to interpolate any Intermediate point (see Fig. 7. 2.8.1). Thus, 
for any parameter <fi 


R-R 


y CL s 


CLH 


^clt’Rclh 


( 7 . 2 . 8 . 1 ) 


and 

= <t>H + 3< £h " ^T )Y CL + {Z 4>h +2< ^T ■' t Avi )Y CL (7-2. 8.2) 


In addition, along the chord line projected to the (r,z) plane a given thickness 
distribution may be applied, according to ISHAFE, using Function FTHICK. 
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7.2.9 Subroutine GDUCT 


Ob.lect Compute Duct Shape 
Options 

I0PT3=1 Straight Annular Duct 
* =2 Read Duct Shape 

=3 Straight -Wall Diffuser 
=4 NACA Curved-Wall Diffuser 

List of Symbols 

(As needed by user . ) 

Theory 

This subroutine is used to prescribe the duct shape r H (z),r T (z), wall bleed 
m H (Z),m T (Z)j and wall temperature T H (z), T t (Z), as required. Since these functions 
are input, the user programmer may write a subroutine for this purpose or read the 
required information according to I0PT3. In addition, the subroutine computes the 
reference radius and normalizes the variables r and T. The variable m is normalized 
in Subroutine FL0WIN when u r is calculated. 

Input /Output 

The user may program any duct shape and wall boundary conditions as required. 

In addition to any input the programmer may write in this subroutine, the programmer 
has available the number of streamwise mesh points JX. The output of this subroutine 
must be (R(l,K,J), 1=1,3; K=l,2, J=1,JL) and Zl. Note that all variables are 
normalized as shown in the sample subprogram described in the Subroutine GDUCT listing 
and that equally spaced spanwise stations are used. The flow chart (Fig. 7-2. 9-1) 
should be followed in programming. 
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FIG. 7. 2. 9.1. -FLOW CHART FOR SUBROUTINE GDUCT. 
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KSSLE,KSSTE, Leading edge and trailing edge index 
SLE,STE , Leading and trailing edge coordinate (dimensionless) 
ZLEH,ZLET , Axial distance hub leading and trailing edge (dimensionless).* 
ZTEH,ZTET , Axial distance tip leading and trailing edge (dimensionless) 




The intersection of the leading and trailing edge of the blade with the hub and 
tip casing are obtained from Subroutine GBLADE. Then the coordinates- of the hub 'and 
tip boundaries are searched until the proper value of streamwise coordinate S for 
the leading edge and trailing edge of the blade are found. The coordinate index ... 
KSSLE is located just upstream of the blade and the coordinate KSSTE is located just 
downstream of the blade. 
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7.2.11 Subroutine S0LVE 


Object Integrate Equations of Motion 
Options 

I0PT4, Print every I0PT4 Station 


List of Symbols 



AAR 

= A /A, 

r, ^ y- 1 2 - 

r*' , 


J4. I 

1 + — 

2(y-l) 

AARS 

= A/A = 




Mx. 

X" 1 


ALPHA 

= a 

L 2 J 

> 

AMACHI 

= Mj. 


> 

c ■ 

= ari(o)/ds 

, V • 

3 

CFSH 

= c /sh e W°)^ 

3 

CFST 

= c, ST * r gn ( 1 ) /zj ( 

3 

CFPH 

= C J<£h = ^ n <0)/q ( 

3 

CFPT 

= C i4>J- C </,n (l)/ ^| 

3 

CPINC 

= C pI =l-(A/A) 2 

3 

CFCOMP 


-n r (o)} 

PRCGEF 

= c p = (P-P(O)) /q, 

3 

QKK 

= q,/(l/Zp T u^ 

3 

Theory 




The equations 

of motion in Section 3* 


the procedures described in Section 3.6. 


Area ratio (dimensionless) 

Critical area ratio (dimensionless) 

Swirl angle (deg) 

Inviscid flow Mach number (dimensionless) 

Constant of integration (dimensionless) 

Streamwise wall friction coefficient hub 
(dimensionless) 

Streamwise wall friction coefficient tip 
(dimensionless ) 

Tangential wall friction, coefficient hub 
(dimensionless) 

Tangential wall friction coefficient tip 
(dimensionless) 

Incompressible pressure coefficient 
(dimensionless) 

Compressible pressure coefficient 
(dimensionless ) 

Local pressure coefficient (dimensionless) 
Average inlet dynamic pressure (dimensionless) 

’ are integrated numerically according to 
Che flow chart is given in Fig. 7.2.11.1. 
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INTEGRATE IN 
STREAMWISE 
DIRECTION 


UPDATE 

SOLUTION 


WRITE 

SOLUTION 



FIG. 7.2.11.1.-FLOW CHART FOR SUBROUTINE S0LVE. 



















7.2.12 Subroutine SPLDC 


Object Compute Constants for Inner Layer Solution 
Options 

None . 

List of Symbols 


UU(1) = u * 

, Velocity (dimensionless) 

UU(2) = H a * 

, Stagnation enthalpy (dimensionless) 

UU(3) = H* 

, Stagnation enthalpy (dimensionless) 

UU(k) = Y + 

, Length (dimensionless) 

UU(5) = W* 

, Mass flow (dimensionless) 

Theory 



This subroutine numerically integrates Eqs. (3.^.5) through (3.*+. 7) out to a 
value of Y m using a third order Range -Kutta method. At the value of Y M> the numerical 
solutions are matched to the analytical solutions given in Section 3*^ and the constants 
of integration are evaluated. 
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7.2.13 Subroutine STRESS 
Object Compute Stress and Heat Flux 
Options 

N0PT2=1 Adiabatic Wall^ 

=2 Nonadiabatic Wall 
N0PT1=1 P + = m + = 0 
=2 P + / 0, m + = 0 
=3 P + =0, mV 0 
N0PT6=1 Centerbody 

=2 No centerbody 

List of Symbols 


AMEH,AMPT 

= 

mn ,mj 

9 

Mass flow bleed hub, tip (dimensionless) 

AMUO 

= 

00 

9 

Free -stream eddy viscosity (dimensionless) 

AMUM 

= 

WfOco 

9 

Free-stream molecular viscosity (dimensionless) 

AWH,AWT 

= 

Q WH’ a WT 

9 

Wall swirl angle, hub, tip (dimensionless) 

DELQ 

= 

A*oo 

9 

Mean displacement thickness (dimensionless) 

KMM 



9 

Index for maximum velocity (dimensionless) 

PFHjPPT 

= • 

Ph ,Pt 

9 

Pressure gradient hub, tip (dimensionless) 

QWH,QWT 

= 

q h 1 Q T 

9 

Wall heat flux hub, tip (dimensionless) 

RH0M 

= 

. P M 

9 

Mean free-stream density (dimensionless) 

SKMH,SKMT 

= 

^•WH ’ ^ HT 

9 

Wall stress hub, tip (dimensionless) 

SWH , SWT 

= 

I(kmh)£(kmt) 

9 

Wall stress matching point (dimensionless) 

TWH,WT 

= 

^WH ,w WT 

9 

Wall temperature (dimensionless) 

t'/M 

= 

U C0 

9 

Maximum free-stream velocity (dimensionless ) 


USTARH, WSTART = U* , U* » Friction velocity hub, tip (dimensionless) 
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Theory 


The inner layer solution is a function of U*(S).- Since the inner layer solution 
does not correct for dll*/dS, U s and i/' are continuously updated to match the inner'and 
outer layer mass flows. 

The normal derivatives, stresses, and heat flux are computed using central 
differences in the core flow. For the inner layer solution, derivatives are obtained 
.analytically from FUPLUS and FHFLUS. 
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7.2.14 Subroutine VARB(PI3,AX,EX,K,I X) 
Object Compute Blade Force Variables 


Options 


IX = 1,2,3 Steps in third order Runge-Kutta method. 


List of Symbols 
A 

AX = dn 3 /d r) 
BX = d\j/ 3 /d77 
IX 
K 

PI3 = n 3 
RH03 = P 3 

A 

TH3 = ® 3 
UPHI3 = U^) 3 
US3 = U S3 
Theory 


, Normal pressure gradient (dimensionless) 

, Normal mass flow gradient (dimensionless) 

, Runge-Kutta index (dimensionless) 

, Streamwise station (dimensionless) 

, Downstream static pressure (dimensionless) 

, Downstream density (dimensionless) 

, Downstream temperature (dimensionless) 

, Downstream swirl velocity (dimensionless) 

, Downstream streamwise velocity (dimensionless) 


When I0FT2=2, the blade force is known, but the loss and exit angle are not 
known. This subroutine calculates the integrands for Runge-Kutta integration of 
Eqs. (3.3.28) and (3*3.29). For a given n 3 , have for iteration i, 


fi, • A t" 


(7.2.14.1) 


A A r 

®3= ®02''l' + 




(7.2.14.2) 
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(7.2.14.3) 


A A A 

p 3 ■ n s /® 5 


s c 


a 2 A f "o i a 2f”: 

U S3~ U S2 L< — ^■J U S3 “ -£~ L“G 


A A 

ft 3 ~ ft 2 

/M r 2 


]- 


A 

P® a 2 


.2 * 

P 3 


A ^S2 ' G 


A A 


2 H 


* 


Ui3 = Ui2 + i-a a 5 a T" 
G B l.P3 U S3 + ?2 4 Is2J 


A(i-I) pa' 2 a 2-11/2 / 1/2 


a<i-I) pa' Z a Z-i 

*3 =M r[u S3 + ^ 3 J 


(7.2.14.4) 


(7.2.14.5) 


(7.2.14.6) 
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7.3 List of External Function Subroutines 


Name Object 


7.3.1 

AMU 

Compute Molecular Viscosity 

7.3.2 

FHPLUS 

Compute Thermodynamic Functions of Inner Layer 

7-3.3 

FUPLUS 

Compute Functions of Inner Layer 

7.3.4 

QWALL 

Compute Wall Heat Flux 

7.3.5 

SPLD 

Functions for Calculating Inner Layer 

7.3.6 

TAWALL 

Compute Adiabatic Wall Temperature 

7.3.7 

USTAR 

Compute Friction Velocity . •••.'' ■< - 

7.3.8 

FTHIK 

Compute Blade Thickness 



11.0 



:• 7.3.3- Function AMU(T ) $.\1 - 

Ob.lect Compute Molecular Viscosity V.T- 

Options •#?. i,-! .v • •, 

None v ;■ 7, w-.: ' 

List of Symbols v - j '■ v , •••*; ,..v. ,- ;f . 

■' >' .• " S- ^ .• v 

AMU Ratio of Molecular Viscosity <■, (dimensionless ) 

-'i V t A*. . .* 

T = 0 , Static Temperature Ratio (dimensionless) 




Theory 

The molecular viscosity is computed according to Sutherland's formula (Ref. 14). 
The working fluid is assumed to be air. Accordingly, ' 


H- ~ 3/2 I + 198.0/ T r 

TT— a 0 " 

r-r 0+ |98.0/Tr 


(7.3.1) 


im 



7-3-2 Function FHPLUS (Y PLUS ,UBL , DUFL , T , DT , SQ, LL ) 
Ob.lect Compute Thermodynamic Functions of Inner Layer 
Options 


N0PT1=1 p + = 0, m + =0 

=2 pV 6, m + = 0 

=3 p + - 0, mVo 


N0PT2=1 

Adiabatic Wall 



=2 

Nonadiabatic Wall 


. ll 

=1 

FHPLUS = Tyy — T a + or 

T + - T + 
'w 


=2 

FHPLUS = dT A /dY + or 

dTVdY* 


=3 

FHPLUS = Q A or 

0 + 


=4 

FHPLUS = dO A /d/ or 

dQ + /dY + 

List of 

Symbols 






DT 

dTVdY + 

DUPL = 

d u + /d y + 

FHPLUS 


UPL = 

u* 

Q 

SQ 

Cx-Om? P w u 

T 

T w -T + 

YPLUS = 

Y* 


, Derivative of temperature (dimensionless) 

} Derivative of velocity (dimensionless) 

, See option list (dimensionless) 

„ Universal velocity dimensionless),.’ , • - • 
, Heat flux parameter (dimensionless) 

, Temperature difference (dimensionless) 

, Universal distance from wall (dimensionless'' 


Theory 


An analytical solution exists for the thermodynamic functions for the inner 
layer as given in Section 3.4. Therefore, the function and its derivatives are known 
analytically according to the appropriate options. 
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7.3-3 Function FUPLUS(UST,YA,LL) 


Object Compute Functions of Inner Layer Solution 
Options 


N0PT1=1 

p* = o,m* = o 

=2 

p + * o ,m + = o 

=3 

p* = 0 ,mV 0 

H 

ii 

FUPLUS = du7dY + 

=2 

FUPLUS = U + 

=3 

FUPLUS = ( yy) 

=4 

FUPLUS = W* 


List of Sv 

mbols 


FUPLUS = 

- 

, See Option List (dimensionless) 

USTA 

U* 

, Friction Velocity (dimensionless) 

YA 

Y 

, Distance from Wall (dimensionless) 

Theory 




An analytical solution for the inner layer U + (Y*) is known from Section 3 A • 
Therefore, the solution and all its derivatives are known analytically. This function 
subroutine calculates the required functions according to the options listed above. 
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7-3.4 Function GWALL(W,U,USTA,R0 w) 

Object Compute Wall Heat Flux . 

Options 

N0PT2=1 Adiabatic Wall * . 

=2 Nonadiabatic Wall 

IF (R(1,3,J) >0.001) N0PT2 = 2 For Outer Walls' 

IF (R(2,3J) >0.001) N0PT2 = 2 For Inner Walls 

List of Symbols 


T = ® 

™ = ® w 
QWALL — ~ Q^f 

U = U 
USTA = U* 
R0W = P 

Theory 


, Static temperature at matching point (dimensionless) 

, Wall static temperature (dimensionless) • 

, Heat flux (dimensionless) = ■ i 

, Magnitude of absolute velocity at matching point (dimensionless) 
, Friction velocity (dimensionless ) m . 

, Density at the wall (dimensionless)' - 


If the wall static temperature is specified, theij the .program assumes that the 
wall heat flux should be calculated. Since the flow conditions are known at the 
matching point, and the friction velocity is known, the enthalpy H*can be calculated 
from Eq. (3-4.25). The heat flux is then given by 


Qw = [®w - ( 


y - 1 


2 2 
M f U 


>] 


P U« 
1 yv 


(7. 3. 4.1) 
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7 . 3.5 Function SPLDtYjU.Hji'jAMjAKjAAjPLjFT) 


Ob.ject Provide Functions for Calculating Inner Layer 
Options 



I =1 

SLPD = 

du + (mVdY + 


=2 

SLPD = 

dH A /dy + 


=3 

SLPD = 

dH + /dY + 


=4 

SLPD = 

— 


=5 

SLED = 

du + (P + )/dY + ‘ - 

List 

of Svmbols 


AA 

= 

A + 

, Van Driest constant (dimensionless) 

AK 

= ' 

K 

, Von Karman constant (dimensionless) ; 

AM 

= 

m + OR 

p + , Mass flow or pressure gradient depending on I (dimensionless) 

H 

= 

H A + OR 

H + ■ , Stagnation enthalpy depending on I (dimensionless) 

PL 

= 

P RL 

, Laminar Prandtl number (dimensionless) 

FT 

= 

P RT , 

, Turbulent Prandtl number (dimensionless) 

SPLD 



, See option list (dimensionless) 

U 

’ = 

U + 

, Universal velocity (dimensionless) 

Y ' 

’ ' = ' 

Y 

,• Universal distance ' (dimensionless ) 


Theory 

This function subroutine provides the. required integrands for integrating 
Eqs. (3.4.5), (3.4.6), and (3.4.7) using a third order Runge-Kutta method. 
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7.3.6 Function TAWALL(T,U,USTA) 
Ob.lect Find Adiabatic Wall Temperature 


Options 


N0PT2=1 Adiabatic Wall 

=2 Nonadiabatic Wall 

IF (R(1,3,J)< 0.001) 
IF (R(2,3,J) <0.001) 


N0PT2 = 1 For Outer Wall 
N0PT2 = l For Inner Wall 


List of Symbols 


T =0 
TAWALL = ® w 
U = U 
USTA = U* 


, Static temperature at matching point (dimensionless) 

, Adiabatic wall static temperature (dimensionless) 

, Magnitude of velocity at matching point (dimensionless) 
, Friction velocity (dimensionless ) 


Theory 


If the -wall static temperature is not specified, then the program assumes that 
the adiabatic wall temperature should be calculated. Since the flow t conditions are 
known at the matching point and the friction velocity is known, the wall temperature 
is obtained from the relation 


®AW 


®+ + U* 2 H*] 


(7. 3.6.1) 


where is determined from Eqs. (3.4.16) or (3.4.20). 
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7.3.7 Function USTAR(KB,YB) 


Object Find Friction Velocity 
Options 

N0PT1=1 P + = 0 ,m + = 0 
=2 pV 0, m + =0 
=3 P +: 0. mVO 

p*> 0 ,IF(p + < 0) p*=0 

rn* > - 0.06; IF (rn*< -0.06) rn* = -0.06 

List of Symbols 


DFE 

= 

dF/du* 

DUPLUS 

= 

duVdv' 

EE 

= 

FVu* 

FE 

= 

F* 

EES 

= 

• e 

ITEM, 



KB 



UPLUS 

= 

u + ' 

USTA 

= 

*(i) 

u 

US TAR 

= 

u* 

YB 

= 

Y 

Theory 




, Derivative of error function (dimensionless) 

, Derivative of velocity (dimensionless) 

, Error (dimensionless) 

, Error function (dimensionless) 

, Error tolerance (dimensionless) 

, Maximum number of iterations permitted (dimensionless) 
, Matching point (dimensionless) 

, Universal velocity at matching point (dimensionless) 

, Guess for friction velocity (dimensionless) 

, Friction velocity (dimensionless) 

, Distance from wall (dimensionless) 


The friction velocity is obtained by matching the inner layer analytical 
solution to the numerical solution at the point KB at a distance YB from the wall. 
This matching requires an iterative procedure since Eqs. (3.4.14), (3*4. 18), or 
(3.4.22) are too complex to solve for u*. The slope of the function U + (Y + ) 4s 
always known analytically so Newton's method is chosen. Let 
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u*= u/u 


( 7 . 3 . 7 . 1 ) 


Y + =, 


n r Pw 


U*Y 


P 


i \ v f an \ 

X N R M r 2 W / p w u * 3 vas /w 


o 



Then an error function is defined as 


F + (U*) = U* (Y*) - u/ U* 


and 


dF* _ N R P w _dl/ u 

dU* (/iw/ff-r) dY + U * 2 


The iteration proceeds by setting 


U* (H) = U* (l) -F(U* 2 )/dF(U* (i) ) /U* (l) 


( 7 . 3 . 7 . 2 ) 


( 7 . 3 . 7 . 3 ) 


( 7 . 3 . 7 . 4 ) 


( 7 . 3 . 7 - 5 ) 


( 7 . 3 . 7 . 6 ) 


( 7 . 3 . 7 . 7 ) 
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Convergence occurs when 



7.3.8 Function FTHIK(Z,IS,L0P) 
Object Compute Blade Thickness Distribution 


Options 


IS=ISHAHS Blade Shape Index, 

L0P=1 FTHIK = t/c 

=2 FTHIK = d(t/c)/d(x/c) 

. fr!. 

List of Symbols 

Z = X/C j Fractional Chordwise Distance (dimensionless) 

Theory 

A parabolic arc thickness distribution of the form 

-b 4 [HH-(-c-) 2 ] (7.3.8. 1) 


is assumed for IS=1. 
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7.4 

List of DIMENSI0N Variables in Blank C0MM0N 

AINFUT(1.J,K) - Y 

, Spanwise location (dimensionless) 

AINRJT(2,J,K) = n 0 

Total pressure (lb/ft abs) 

ainput( 3 ,j,k) = n 

Q 

, Static pressure (lb/ft^ abs) 

ainput(4,j,k) = a 

, Svirl angle (deg to axis) 

AIKPUT(5,J,K) = © Q 

, Total temperature (deg R) 

J = 1 

, Upstream of blade row (dimensionless) 

1=2 

, Downstream of blade row (dimensionless) 

K = 1,KLL , Number of spanwise stations 

BINPUT(1,J,K) = Y 

, Spanwise location (dimensionless) 

BINRJT(2,J,K) = n 0 

2 

, Total pressure (lb/ft .abs) 

binfut( 3 ,j,k) = n 

Q 

, Static pressure (lb/ft abs) ' 

binfut(4,j,k) = a 

, Swirl angle (deg to axis’! 

BINPUT(5,J,K) = ® 0 

, Total temperature (deg R) 

J = 1 

, Inlet flow 

OJ 

II 

*"3 

, Exit flow 

K = 1,KLL 

, Number of spanwise stations 

C0NST(l,L) C, (m + ) 

, Velocity (dimensionless) 

C0NST(2,L) C 2 (m + ) 

N 

Adiabatic stagnation enthalpy (dimensionless) 

C0NST(3,L) C 3 (m + ) 

, Nonadiabatic stagnation enthalpy (dimensionless) 

C0NST(4,L) C 4 (m + ) 

Stream function (dimensionless) 

C0NST(5,L) C 5 (p + ) 

, Velocity (dimensionless) 

C0NST(6,L) C 6 (p + ) 

, Stream function (dimensionless) 
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CNST(l) = 

l =: 

DF(J,K) = 
DHF(1,K) = 
DHF(2,K) - 
DFF(l,K) = 
DPF(2,K) > 
K : 

F(2,I,K) ■ 
F(3 jl s K) i 
F(M,K) ■ 
F(5 jIj'K) ■ 
F(6,I,K) ■ 

f( 7 jl » k) 

I ! 


K 

FF(1,I,K) 

FF(2,I,K) 

FF(3,I,K) 

ff(4,i,k) 

FP(5,I,K) 


C0NST(J,L) • Interpolate values between L, ‘L+l 


dF (j, 2, K j/d'S t Streamwise derivatives (dimensionless ) 


(fy/dS)" 

(an/as) H 

(a^/as) P 

(an/ds) P 

1,KL 

$ 

u s 

U<£ 

n 


i 

© 

p 


, Hbmogeneous solutidn 'for stream function (dimensionless) ' 
, Homogeneous solution for static pressure (dimensionless) 

, Particular solution for stream function (dimensionless) 

, Particular' solution for static pressure (dimensionless) 

", Number of streamlines 

,'Stfeam' function' (dimensionless) ' ' 

, StreSirwise'' velocity ''(dimensionless) 

-, : Swijl- velocity., (dimensionless),- f v ~.y- r 

, Static -pressure (dimensionless) « 

, Entropy (dimensionless)^.- : • 

/Static . temperature, (dimenionless.) - - 

, Density (dimensionless 


1 


, Inlet conditions - 


= 2 


= 3 


= 1, KL; 

A 

= M ' 

A 

= n 

A~ 

= © 

" ®6 - 
= Ao 


9 S -S \ ... r ■ ./) . ^ 

jS^S + dS'-' 

Number- of . streamlines-, » .7 

VlnviscidiMach- number (dimensionless) . 

^ Inviscid static pr essure-ddimansionless ) 

, Inviscid static: temperature^ (dimensionless) , 

, Inviscid total-.vtemperature - (dimensionless) 

, Inviscid totalJpressure (dimensionless) r 
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ff(6,i,k) = P 

FP(7,I,K) = U s 
FF(8,I,K) - U<£ 

fp(9.»i,k) - 

FF(10,I,K) = U B 
FF(11,I,K) = 3 
FF(12,I,K) = $ 

FF(13,I,K) = I 
FF(l4,I,K) = U 
FF(15,1,K) = z B 

FF(15,2,K) = AI b 
~FF( l6,I,K) = $ 
FF(17,1,K) = C L 
FF(17,2,K) = C D 
1=1 
= 2 

K = 1, KL 
fg(IjK) = a 
FG(2,K) = n 0 
FG(3,K) = ©0 
FG(4,K) = M 
FG(5,K)\= p 0 
FG(6,K) = U 

K = 1, KL 

/ 


, Inviscid density (dimensionless) > f •• 

, Inviscid streamwise velocity (dimensionless) 

, Absolute swirl velocity (dimensionless) c . • - 

, Relative swirl velocity (dimensionless) , 

, Blade velocity ( dimens ipnless) . ....... 

, Absolute angle to axis (deg) . 

, Relative angle to axis (deg) ... . _ . • . .... 

j Inviscid flow entropy (dimensionless) 

, Magnitude of relative inviscid flow velocity (dimensionless) 
j Loss coefficient (dimensionless) 

, Blade entropy rise (dimensionless) K 

, Stream function (dimensionless) 

' , Lift coefficient(dimensioiiless) 

, Drag coefficient (dimensionless) 

, Upstream of blade row 
, Downstream of blade row ■ 

, Number of streamline 
, Inlet swirl angle (deg) 

, Inlet stagnation pressure (dimensionless ) 

, Inlet stagnation temperature (dimensionless). 

»• Inlet Mach' number (dimensionless) - . '• 

> Inlet stagnation density (dimensionless) 

Inlet magnitude of velocity (dimensionless) " 

Number of .streamlines * ■ ~ 
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F0RC(1,K) 

= H s 

F0RC(2,K) 

= H<£ 

F0RC(3,K) 

= a s 

F0RC(4,K). 

= 

F9*RC(5,K) 

= % 

K 

= 1,KL 

G(1,K) = 

/ M-r 

G(2,K) = 

£sn 

G(3,K) = 

£</>n 

G(4,K) ■ = 

dU s /dft 

G(5,K) = 


G(6,K) = 


G(7,K) = 

d ^ /d71 

G(8,K) = 

U^)-Ub 

g(9>k) = 

u 

G(10,K) = 

Esn 

G(U,K) = 


o(ia,K) = 

d@/d% 

G(14,K) = 

di/dn 

G(15,K) = 

Q 

G(16,K) = 

dQ/dfl 

QU,K) = 

R 

Q(2,K) = 

Z 


, Streamwise force/area (dimensionless) 

, Swirl force/area (dimensionless) 

, Streamwise force /span (dimensionless) 

, Swirl force/span (dimensionless) 

, Blade dissipation /area (dimensionless) 

, Number of streamline (dimensionless) 

, Ratio of turbulent to reference molecular viscosity (dimensionless) 
, Streamwise stress (dimensionless) 

, Tangential stress (dimensionless) 

, Normal derivative of streamwise velocity (dimensionless) 

, Normal derivative of' swirl velocity’ (dimensionless) 

, Normal derivative of streamwise stress (dimensionless ) 

, Normal derivative ' of tangential stress (dimensionless) 

, Relative velocity (dimensionless) 

, Magnitude of relative velocity (dimensionless) 

, Streamwise rate of strain (dimensionless) 

> Tangential rate of strain (dimensionless) 

, Normal derivative of static temperature (dimensionless) 

, Normal derivative of entropy (dimensionless) 

, Heat flux (dimensionless) 

, Normal derivative of heat flux (dimensionless) 

, Radius (dimensionless) 
j Axial Distance (dimensionless) 
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r 


Q(3,K) = dR/d% 
q(4,k) = aR/as 
q(5,k) = a 2 R/a^/as 
q(6,k) = v 
Q(7,K) = . av/dtf 
q(8,k) = av/as 
q( 9 ,k) = a 2 v/a?j/as 

Q(10,K) = Y 
Q(11,K) = Y/Y t 
Q(12,K) = A 
Q(13,K) = G 
q(i4,k) = aG/a?z 
q(i5,k) - aG/as 
q(i6,k) = aya?7 
Q(17,K) = a 2 7 ]/d?l 2 
Q(18,K) = 71 
Q(19,K) = 77 
K = 1,KL 

R ( I j 1 j J ) = R T ( z j ) 
R(2,l,j) - Rh(Zj) 
R(l,2,j) = m T (z T ) 
R(2 ,2 , j) = m H (Zj ) 
R(1,3,J) = ® H (Zj) 

R(2,3,J) = ® T (Zj) 


, Normal derivative of radius (dimensionless) 

, Streamwise derivative of radius (dimensionless) 

, Second derivative of radius (dimensionless) 

> Metric scale coefficient (dimensionless) 

) Curvature of potential line (dimensionless) 

> Curvature of streamline (dimensionless) 

> Second derivative of metric scale coefficient (dimensionless ) 

> Physical distance across duct (dimensionless) 

, Fractional distance across duct (dimensionless ) 

, Area between ajacent streamlines (dimensionless) 

, Gap between blade surfaces (dimensionless) 

, Normal derivative of blade surface (dimensionless ) 

, Streamwise derivative of blade surface (dimensionless ) 

, "Transform of normal coordinate (dimensionless) 

, Second derivative (dimensionless) 

, Normal coordinate (dimensionless) 

, Transformed normal coordinate (dimensionless) 

5 Number of streamlines (dimensionless) • 

/ 

» Radius of hub (dimensionless) 

, Radius of tip (dimensionless) 

, Mass flow of tip bleed (dimensionless) 

, Mass flow of hub bleed (dimensionless) 

, Wall temperature of tip (dimensionless) 

Wall temperature of hub (dimensionless) 
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R0(1,K) = 

R 

, Radius (dimensionless) 

R0(2,K) = 

dR/dz 

, Derivative of radius (dimensionless) 

R0(3,K) = 

d 2 R/dz 2 

, Second derivative of radius (dimensionless) . 

R0(4,K) = 

Z 

, Axial distance (dimensionless) 

R0(5,K) = 

V 

, Metric. scale coefficient (dimensionless) 

R0(6,K) = 

dY/dS 

: , Derivative of metric scale coefficient (dimensionless) 

R0(7,K) = 

Yr 

, Distance .across duct (dimensionless) 

R0(8,K) = 

s 

.j Streamwise coordinate (dimensionless) , 

R0(9,K) = 

m 

, Mass flow bleed (dimensionless) 

R0(io>k)= 

©w 

j Wall temperature (dimensionless) . 

R0S(l ) = R0(l ,K) 

. , . . - .... , ■ 


Note: 3 arrays defined where 0 takes on the value,. .H,. M, T 

: H , Hub Wall ' ; - - 

ivj" -Me’an line' " ' • ' r • 


T ,Tip wall 



7.5 List of Variables in Blank C0MM0N 


ACHI 

c 

X 

, Clauser constant (0.016) ' 

ACHAPPA 

= 

K 

, von Karman constant (0.4l) ' 

ALP1 

= 

Q l 

, Hub swirl angle (deg) 

ALPHS 

= 

a s 

, Stagger angle to axis (deg) 

ALPSH 

= 

a H 

, Hub stagger angle to axis (deg) 

ALPSM 

= 

a M 

, Midpoint stagger angle to axis (deg) 

ALPST 

= 

“t 

, Tip stagger angle to axis (deg) 

AMACH1 

= 

M, 

, Average inlet Mach number (dimensionless) 

AMACHR 

= 

Mr 

, Reference Mach number (dimensionless) 

AMS1 

= 

Ms. 

* 

» Hub Mach number for input (dimensionless ) 

AMPLUS 

= 

m* 

, Mass flow bleed parameter (dimensionless) 

ANH,ANT 

= 

n H. n T 

, Power law profile (dimensionless) 

APLUS 

= 

A + 

, van Driest constant (26.0) 

AFRES1 

= 

P| 

, Average inlet static pressure (dimensionless) 

AREAR 

= 

Or 

p 

, Reference area (ft ) 

CHORD 

= 

B 

, Strut chord local (dimensionless ) 

CHORDH 

= 

b H 

, Strut chord Hub (dimensionless) 

CHORDM 

= 

b M 

5 Strut chord midpoint (dimensionless) 

CHORDT 

= 

&T 

> Strut chord tip (dimensionless ) 

CPR 

= 

c pr 

5 Specific heat at constant pressure (5997.0 ft /sec^/deg R) 

CPV 

= 

o 

T? 

< 

> Specific heat at constant volume (3283.0 ft 2 /sec 2 /deg R) 

DDS 



> Mesh distortion parameter (1. < DDS < 10. ) 

DETA 

= 

Atj 

, Step size in Normal coordinate (dimensionless) 
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dmplus 

= ■ 

Am* 

, Step size in asymptotic constant table (m + ) (dimensionless) 

dplusm 

= 

AP + 

, Step size in asymptotic constant table (p*) (dimensionless) 

BSH,DST 

= 

A* A* 

, Displacement thickness hub, tip (dimensionless) 

DS 

= 

AS 

, Streamwise step size between stations (dimensionless ) 

DSS 

= 

dS 

, Streamwise step size (dimensionless) 

dytjpi 

= 

di 

, Average inlet dynamic pressure (dimensionless) 

EP 

= 


, 2.7182818 

gamma 

= 

7 

> : r r 

, Ratio of specific heats (1.4) 

GAP 

= 

G 

Gap between blades (dimensionless) 

gasr 

= 

<R 

, Gas .constant (1714.0 ft /sec / deg R) 

GMR1 

= 

(/“ 1) Mf 


GMR2 

= 

yUf 


gravr 

= 

g 

***■.• • p 

, Gravitational constants (32.2 ft/sec ) 

ishape 



, Blade shape parameter (dimensionless) 

JL 



, Number of streamwise stations (dimensionless) 

KDS 



, Number of streamwise steps/station (dimensionless) 

KL 



f. •• • ; . - 

, Number of streamlines (dimensionless) 

KLL 



, Number of input streamlines (dimensionless) 

KLMH 



t ... 

, Hub matching point - inner layer (dimensionless) 

KMT 



, Tip matching point - inner radius (dimensionless) 

LM 



, Size of table of constants for inner layer (dimensionless) 

lmm 



, Midpoint of table of constants for inner layer (dimensionless) 

NB 



t i 

, Number of struts (dimensionless) 
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FHIC 

= 

4>c 

, Blade camber (deg) 

, * 4 

PHICH 

= 

$CH 

, Blade camber - hub (deg) ’ 

PHICM 

= 

■ $CM 

, Blade camber - midpoint (deg) 

PHICT 

= 

$CT 

, Blade camber - tip (deg) 

PPLUS 

= 

: P + 

, Pressure gradient parameter (dimensionless) 

n 

= 

n 

,3-1415926 

FEESSR 

= 

Pr 

, Reference static pressure (PSP) 

PRL 

= 

Prl 

, Prandtl Number laminar 0.70 

PRT 

= 

P RT 

, Prandtl number turbulent 0.72 

RADR 

= 

r r 

, Reference radius (ft) 

REY 

= 

' N r 

, Reynolds number (dimensionless) 

RCLH 

= 

^CLH 

,Hub radius of blade centerline (dimensionless) 

Rem 

= 

r clm 

, Midpoint radius of blade center (dimensionless) 

RCLT 

= 

^CLT 

,Tip radius of blade centerline (dimensionless) 

RHOR 

= 

Pr 

, Reference density (slugs/ft^) 

SL 

= 

$L 

, Length of duct in streamline coordinates (dimensionless; 

SNDR 

= 

Cr 

, Reference speed of sound 1116.0 (ft/sec) 

SOLD 

= 

cr 

, Solidity (dimensiohless ) 

TEMPR 

= 

Tr 

, Reference temperature (deg- RANKIN) 

THICK 

= 

t 

, Local blade thickness (dimensionless) 

THICKH 

= 

♦h 

, Blade thickness hub (dimensionless) 

THICKM 

= 


Blade thickness midpoint (dimensionless) 

THICK 

= 

, -t T 

, Blade thickness tip (dimensionless) 
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THIKM 

= 

t/B 

, Maximum thickness/chord (dimensionless) 

TI 

= 


, (017^5329 radians /deg) n 

USR 

— . i 

. u r 

, Reference radius (dimensionless) 

VI SCR 

= ' 

H-r 

, Reference molecular viscosity (0. 370xlO~^) 

WPLO 

= 

Cl) 

, Weight flow (lb/sec) 

YPLUSM 

' = 

< 

, Matching point for table of asymptotic constants (dimensionless) 

ZCL 

.= ' 

Z CL 

, Axial distance to blade centerline (dimensionless) 

ZLE 

= 

Z LE 

, Blade leading edge (dimensionless) 

ZTE 


*TE 

, Blade trailing edge (dimensionless) 

0 

ZI 

= 

z. 

, Duct axial length (dimensionless) 



.7-6 Description of Mesh Geometry 

Since the numerical solution of turbulent boundary layer equations is required, 
the mesh points for the normal coordinate are not distributed uniformly across, the 
duct, but crowded near the duct walls. This nonuniform distribution of mesh points 
is accomplished with the r)( D ) transformation given by Eq. (7. 2. 4. 5) in Subroutine 
C00R. Thus the actual calculation is done, in an ( 77 ( n ),S ) coordinate system. 

The ( t], S ) coordinate system is shown in Fig. (7.6.1). It consists of KL 
uniformly distributed mesh points on the normal coordinate and JL uniformly distributed 
mesh points in the streamwise direction. Since the explicit numerical integration 
method used in this report* is conditionally stable, the coordinate in the streamwise 
direction usually requires a finer mesh 'then specified by JL. Thus the JL mesh points 
can further divide by KDS . Therefore ' ' 

(7.6.1) 

) 


and 

S =AS (J-I) + dS(KS-l) J= I , JL ; KS = l,KDS (7.6.2) 


17 = A 17 (K-l) K = I , KL 


The transformation 17 (n) involves a parameter C which can be related to the 
distortion of mesh points by 


n DPS 
I +DDS 


(7.6.3) 


where 


DOS 



(7.6.4) 


Finally, it is noted that matching the inner layer solution with the numerical 
solution occurs at KMH and KMT. 
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8.0 APFENDIX B: SAMPLE TEST CASES 


A number of test cases are shown in this section to illustrate the various 
options available in the computer program and to test the accuracy of the calcula- 
tion by comparison with experimental data where possible. Included with the write- 
up is an input data sheet used to set up the- problem. 



8.1 Boundary Layer in Entrance Region of Circular Pipe 

This case illustrates use of the computer program to obtain the development of 
the boundary layer in the inlet region of a circular pipe where the streamwise 
pressure gradient is small and the boundary layer is thin. The calculated wall 
friction coefficients as. a function of Reynolds number based on momentum thickness 
for Mq-, = 0.1 with different values for von Kar man constant k are plotted in 
Fig. 8.1.1 and compared with the wall' friction laws obtained from Ref. l4. Calculated 
wall friction coefficients and recovery factors as a function of Reynolds number 
based on momentum thickness for a range of Mach numbers are plotted in Fig. 8.1.2. 

The computed wall friction coefficients are compared with the correlations of 
Ref. 20. 





i v.-.; 
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RECOVERY FACTOR^ SKIN FRICTION COEFF.ICIENTxlO 


PRESENT COMPUTATIONS 
K = 0.40 

V ' Prt = 0.9 . 

M =c , 

0 ' o . l • 
o 2.0 

□ 4.0 




FIG. 8.1.2r-VARIATION OF SKIN FRICTION COEFFICIENT AND • 
RECOVERY FACTOR WITH REYNOLDS NUMBER. 
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8.2 Conical Pipe Diffuser 


The computed effect of compressibility on the performance of a subsonic conical 
pipe diffuser is shown in Fig. 8.2.1 and compared with the experimental data from 
Ref. 21. The solutions are in good agreement with the experimental data showing the 
increase in pressure recovery with increasing Mach number as well as a drop in 
pressure at the diffuser throat. This trend is also indicated by the data in Ref. 5. 
Shown also on Fig. 8.2.1 is the incompressible solution for this diffuser obtained 
from the numerical procedure described in Ref. 6. These results follow the trend in 
pressure recovery with Mach numbers; however, the trend in the predicted separation 
point with Mach number is not consistent with that for the present compressible 
solutions. This may be partially explained by the fact that the compressible solu- 
tion tends to calculate slightly higher wall friction coefficients than measured 
experimentally (see Fig. 8.1.1). 



PRESSURE COEFFICIENT ^ - P. )/(P 



FIG. 8.2. 1.-EFFECT OF INLET MACH NUMBER ON PRESSURE COEFFICIENT 
IN A 10° CONICAL DIFFUSER. 
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8.3 Annular Diffuser 


The computed pressure coefficients for a curved-wall annular diffuser are 
presented in Fig. 8.3.1 for two different inlet Mach numbers ( = 0.47 and M, = 0.75) 

and compared with the data of Ref. 22 . The predicted pressure recoveries are in good 
agreement with the experimental data and the trend of increasing pressure recovery 
with increasing inlet Mach number is again demonstrated, it should be noted that 
the computed wall pressure coefficients presented in Fig. 8.3.1 are not those obtained 
directly from the computer output but were corrected to allow direct comparison with 
the experimental data which is based on a dynamic pressure obtained from a free- 
stream Mach number rather than a weight -flow -averaged Mach number. 



PRESSURE COEFFICIENT (P 




FIG. 8.3.1 - EFFECT OF INLET MACH . NUMBER ON PRESSURE 
COEFFICIENT FOR AN ANNULAR DIFFUSER. 
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8.4 Curved-Wall Diffuser With Struts 


■ % 

The computed pressure recovery coefficients for the NASA curved -wall diffuser 

( M,= 0.47) described in Ref. 22 and Section 8.3 were recalculated with the. addition 
of two streamline struts placed in the diffuser as .shown in Fig. 8.4;1. These 
streamline -struts were NASA 5 digit series uhcambered airfoil sections with a one 
foot chord and 20 percent thickness -to -chord ratio. The primary effect of these 
struts is to lower the static pressure coefficient in the region of the strut due to 
the increased flow Mach number. There has also been a slight loss in pressure 
coefficient due to the blade loss. It should be noted that the solution does not 
account for the slight additional flow blockage due to the strut boundary layer 
growth. In addition, the increased losses due to the interaction of the strut 
boundary layer with the annulus wall boundary layer are not taken into account . 






.f 


l44 


\ 



1.2 

0 1.0 2.0 3.0 4.0 

DISTANCE FROM DIFFUSER THROAT, X/D 


FIG. 8.4.1. - COMPUTED STATIC PRESSURE DISTRIBUTION IN AN ANNULAR 
DIFFUSER WITH STREAMLINE STRUTS. 
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'8.5 Curved -Wall Annular Duct With Bleed 

In Section 8.3, the flow through an annular diffuser at M, = 0.47 was found to 
separate at a distance X/D = 3*5 (see Fig. 8.3.1). This separation occurred on the 
tip. wall. Since blowing tends to promote separation and suction tends to delay 
separation, a solution was obtained for blowing on the hub wall and suction on the 
tip wall with bleed rates on each wall of 0.02 lb/sec/ft^. Since there is 
essentially no net change in mass flow with and without bleed for this diffuser, a 
comparison of the velocity distributions can be made without compensating for 
differences in mass flow. This comparison, shown in Fig. 8.5.1, indicates that the 
hub boundary layer is blown out from the wall and the tip boundary layer is sucked 
in towards the wall. In addition, the computations predicted that separation on the 
tip wall was delayed to X/D= 3*7. 



DIMENSIONLESS STREAMWISE VELOCITY, 


0»02 

m_j. = 0*02 



DIMENSIONLESS DISTANCE, Y/Y ^ 


FIG. 8.5.1. -VELOCITY DISTRIBUTIONS. IN A CURVED-WALL DIFFUSER 
WITH BLOWING AND SUCTION. 
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